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ABSTRACT 


The  problem  of  narrow  band  interference  while  transmitting  broad-band  signals 
like  Direct  Sequence  Spread  Spectrum  is  a  common  source  of  problems  in  Electronic 
Warfare.  This  can  occur  either  due  to  intention£il  jamming  or  due  to  unavoidable 
signal  sources  present  in  the  vicinity  of  the  receiver.  Lack  of  improper  information 
on  these  narrow  band  interferers  makes  it  difficult  to  cancel  them. 

In  this  thesis  the  above  problem  is  addressed  by  using  an  adaptive  notch  filtering 
technique.  Before  adopting  such  a  technique  other  methods  like  the  Least  Square 
Estimator  and  the  Maximum  Likelihood  Estimator  were  explored.  However  the  Kwan 
and  Martin  adaptive  notch  filter  structure  was  found  both  relevant  and  suitable  for 
the  problem  of  interest.  The  Kwan  and  Martin  method  has  the  difficulty  of  increasing 
hardware  complexity  with  number  of  notches.  This  makes  it  difficult  to  implement 
in  real  time.  A  new  algorithm  was  developed  for  the  purpose  of  implementing  the 
structure  in  real-time.  This  new  algorithm  offers  the  same  performance  at  reduced 
hardware  complexity.  This  algorithm  was  simulated  and  the  results  were  presented.  A 
hardware  feasibility  is  discussed  by  proposing  a  simple  structure  based  upon  existing 
commercial  signal  processing  chips. 
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I.  INTRODUCTION 


A.  MOTIVATION 

The  problem  of  suppressing  narrow-band  noise  from  a  broad-band  spectrum  is 
of  considerable  importance  in  the  areas  of  Electronic  Warfare  and  Anti  Submarine 
Warfare  scenarios.  In  this  thesis  a  workable  solution  is  proposed  for  the  problem  cre¬ 
ated  by  narrow-band  interference.  This  solution  will  come  in  the  form  of  an  adaptive 
digital  notch  filUr.  However  before  the  details  of  the  solution  to  the  problem  are 
discussed,  a  brief  description  of  the  problem  is  given 

1.  Electronic  'Warfare(EW) 

A  typical  EW  scenario  is  given  in  Fig  1.1.  In  this  scenario  a  transmitter 
transmits  information  over  a  long  distance.  As  an  Electronic  Counter  Measure(ECM) 
this  information  is  coded  and  ifansfnitt^  as  a  Direct  Sequence  Spread  Spectrum 
(DSSS)  signal.  Inherently  this  signal  is  a  broad-band  signal.  However  at  the  £W 
receiver  there  is  often  narrow-band  interference  from  such  things  as  push  to  talk 
sysieins(FTS)  that  swamp  out  the  received  DSSS  signal.  Sonretimes  for  reasons  of 
signal  security  PTS  operating  frequencies  are  varying  with  time.  Under  condilioits 
such  as  these  the  EW  receiver  cannot  function  eflecUvely.  For  proper  fnactioning  of 
the  EVV  receiver  we  must  enhance  the  received  signal  by  selectively  and  adaptively 
suppressing  these  narrow-band  signals. 

2.  Assumptions 

This  thesis  addresses  the  problem  mising  at  the  tactical  data  ctunniunica- 
tion  link.  In  this  EW  scenario  as  dipicled  in  Fig  1.1  we  are  attempting  to  perform 
signal  analysis  on  a  DS^  sigiud  using  a  £W  reciever.  This  EW  receiver  is  not  the  dcs- 
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EW  SCENARIO 


Inter-  Inter-  Inter¬ 
ference  ference  ference 


ignated  receiver  and  hence  does  not  have  the  spreading  code  available.  For  this  reason 
the  narrow-band  interference  severely  limits  the  ability  to  determine  the  charaterstics 
of  the  recieved  signal  such  as  carrier  frequency,  chip-frequency  etc. 

A  typical  DSSS  system  is  the  Tactical  Information  Distribution  Systems 
(JTIDS)[Ref.  7].  Data  on  such  a  link  is  generally  in  the  form  of  digital  messages. 
A  typical  L-band  (960-1215  Mhz)  JTIDS  allows  transfer  of  digital  data  between  any 
properly  equipped  users  within  line-of-sight.  The  typical  RF  band-width  required  is 
around  lOMhz. 

Assuming  a  Superheterodyne  EW-receiver,  the  band-width  of  the  Interme¬ 
diate  Frequency  (IF)  need  be  only  lOMhz.  This  enables  us  to  digitize  the  signal  at 
the  base-band  with  a  (20-50)Mhz  Analog  to  Digital(A/D)  converter.  It  is  assumed 
that  this  signal  is  digitized  and  converted  into  a  floating  point  number.  In  this  thesis 
it  is  assumed  that  the  floating  point  arithmetic  is  of  sufficient  precision  that  rounding 
errors  are  not  a  significant  problem. 

B.  PROBLEM  DEFINITION 

The  signal  described  in  section  (I.A.3)  can  be  thought  of  as  a  broad-band  signal 
with  additive  narrow-band  signals  occurring  at  arbitrary  frequencies  and  times.  These 
narrow-band  signals  can  be  modeled  as  either  pure  sinusoidal  signals,  as  outputs  of  a 
narrow-band  filter  excited  by  white  noise,  or  as  any  band-limited  narrow-band  signal. 

We  shall  assume  that  the  signal  y*  is  received  by  the  receiver  and  is  defined  as 

3 

yk-m-^Y^x\  +  -fk  (1.1) 

fat 

where  x*,  the  interferer,  is  either  a  pure  sinusoid  or  a  narrow-bond  signal  whose 
frequency  (i.-  center  frequency)  can  vary  slowly  with  time.  The  terra  lu*  is  the 
broad-band  DSSS  signal  of  Interest  and  71.  is  White  Gaussian  Noise(WGN)  with  zero 
mean  and  variance  or*  A^(0, 0'')  The  goal  is  to  remow  only  the  narrow-band  interferers 


and  obtain  the  output  {wk  +  7*)  for  further  processing.  A  typical  spectrum  for  yk  is 
given  in  Fig  1.2. 

In  this  thesis,  the  design  for  an  adaptive  digital  filter  that  will  achieve  the  goal 
of  removing  narrow-band  interferers  from  a  broad-band  signal  of  interest  is  developed. 
This  thesis  is  organized  into  five  chapters,  chapter  I  is  the  introduction  The  chapter 
II  concentrates  on  building  necessary  background  on  such  things  as  the  Least  Square 
Estimator  and  Maximum- Likelihood  Estimator.  Chapter  III  briefly  explains  existing 
methods  of  adaptive  filtering  with  emphasis  on  adaptive  notch  filtering  techniques 
.  The  new  algorithm  (Ref.  26)  developed  for  the  purpose  of  addressing  the  current 
problem  of  interest  is  also  explained.  Chapter  IV  is  totally  devoted  to  modelling 
various  signals  required  for  testing.  Simulation  results  are  discussed  and  presented. 
Chapter  V  provides  a  brief  survey  of  existing  hardware  components  and  a  look  at 
possible  hardware  structures  and  their  limitations.  Chapter  VI  is  conclusions. 


gnal  Spectrum 


Fig 


II.  ESTIMATION  TECHNIQUES 


A.  DIGITAL  FILTERS 

Digital  Filters  can  be  characterized  by  their  impulse  response,  their  trans¬ 
fer  function,  cr  by  a  difference  equation  [Ref.  13]  .  A  typical  Inhnite-Impulse- 
Respon3e(IIR)  lilter  represented  in  difference  equation  form  is  given  by 

Xk  —  dlXk-l  +  +  <i3Xk~3  +  biUk  +  b2Uk-l  (2.1) 

while  an  Finite  Impulse  Response(FIR)  filter  is  given  by 

Xk  —  biUk  +  b^Uk-i  +  biUk~3  (2.2) 

where  Uk  is  the  input  to  the  iiiter  and  Xk  is  the  output  of  the  filter.  Xk  and  Uk  are 
the  sampled  values  of  the  continuous  fun-lions  x(t)  and  u(t).  By  choosing  as  the 
sampling  time  we  notionauiy  write  the  discrete '  ‘gnal  as 

Xk^x{k)^x{kT,)  (2.3) 

Let  A'(D)  be  the  spectrum  of  the  signal  Xk^  where  H  is  the  angular  frequency 
component  present  in  the  actual  signal.  For  analysis  purposes  consider  normalized 
frequency  given  as  <n  s  (ITf  With  this  definition  we  do  not  ha.e  to  use  the  actual  fre¬ 
quencies  but  only  the  normalized  frequencies.  Asiuming  that  there  is  no  aliasing  while 
sampling  the  maximum  normalized  frequency  is  0.5  cycle/sample  or  jt  rad/sample. 

1.  Some  Hardwcurc  Schemes 

Fig  2.1  represents  a  conventional  Direct  Form  II  t -presentation  of  a  digital 
filter  [Ref.  18).  However  for  high-speed  throughput  it  is  often  necessarj*  to  have  a 
pipelined  architecture.  FIR  filters  are  highly  amenable  for  pipelining  thus  giving  a 
high  throughput  rate  for  computations. 
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Direct  Form  11 


2.  Modelling  Hardware  Multiplier  and  Adder 

Multiplication  of  two  variables  say  6i  and  Uk  in  hardware  can  be  modelled 
as 

b\*Uk  =  bi*  UkDi  =  biUkDi  (2.4) 

where  i  is  the  hardware  multiplier  and  *  is  the  ideal  multiplier  and  Di  is  the  propa¬ 
gation  delay  associated  with  the  multiplication.  Similarly  a  hardware  adder  is  desig¬ 
nated  as  +  while  the  ideal  adder  is  given  as  +  and  the  corresponding  delay  associated 
with  the  adder  is  Dj.  Then  an  addition  of  two  variables  tt*  and  Uk~i  can  be  written 
as 

Ufc+ujb-i  =  (uk  +  tifc-j)D2  (2.5) 

Each  of  the  previous  devices  can  be  incorporated  into  a  pipeline  digital  filter  structure 
by  following  it  with  a  synchronizing  register.  The  clock  period  must  then  be  less  than 
i)  4*  +  tr  where  t,  is  the  register  setup  time  and  U  the  register  propagation  delay. 

3.  Pipelining  an  FIR  filter 

In  this  section  wc  wish  to  realize 

xjt  =  ♦  Ufc  +  •  ut-i  -f  63  *  Ui-j  (2.6) 

It  is  acceptable  to  have  an  overall  delay  ,  However  it  is  desirable  to  minimize  N: 

i*  sa  •ut  +  6a*  Ui-.j  +  fej*  u*»3}  (2.7) 

First  associate  one  D  with  each  multiplier  and  form  real  multipliers  by  the  notation 
D*  =  *.  Then  we  gel 

ik  —  (6ijD  ♦  ua  -f  63/?  •  uA*  1  +  b^D  *  ua-a}  (2-S) 

Now  form  the  real  multiplier 

X*  =i  {6t*ua  +  biiuk  -  1  +  63*ui,.a}  (2.9) 
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Now  we  use  one  additional  D  to  associate  with  the  adding  of  b2*‘>^k-i  and  bz*Uk-2. 
(Note:  It  is  important  to  choose  the  most  delayed  values  first  in  order  to  minimize 

N) 

Xk  =  {Dh*Uk  +  Dib^iuk-i  +  bz*Uk-2)}  (2.10) 

Now  form  the  reed  adder 

Xk  =  {Dbxiuk  +  (b2*Uk,k+hink.2)}  (2.11) 

Now  we  use  an  additional  delay  D  to  associate  with  the  final  adder  to  get  a  real 
adder; 

Xk  =  D'^~^  {Dbi*Uk+(h2*Uk^Ab3*Uk.2)}  (2.12) 

Finally  we  choose  the  system  delay  such  that  D  =  and  replace  Uk-2  = 
z~^Uk  and  Uk-\  —  thus  we  get 

Xk  s=  |j"'6iiu*4-(s"*6jiuj84-s"*63iui)j  (2.13) 

Now  factor  out  of  the  expression  to  obtain 

Xk  =  s''^‘*^{bi*xikHb2iuk^z~%*Uk)}  (2.14) 

Note  that  A' s  2  is  the  minimum  possible  value  of  N  for  a  causal  filter.  Choose 
A^  ~  2  then 

ss  {6,iu|,+(62iUk4-C’’fc3iu*))  (2.15) 

Equation  (2.15)  is  very  convenient  for  pipelining  and  is  given  one  to  one  in  Fig  2.2. 

The  above  procedure  is  general  for  any  FIR  filler.  Pipelining  an  FIR  filler 
results  in  maximizing  the  sampling  or  throughput  frequency  at  the  expense  of  a  delay 
or  latency  in  the  availability  of  the  output. 
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4.  Pipelining  HR  filters 

Pipelining  an  HR  filter  is  considerably  more  difficult  than  an  FIR  filter. 
This  is  mainly  because  the  delay  in  the  availability  of  the  output  makes  it  impossible 
to  feed  back  output  values  with  short  delay  as  is  required  for  the  straightforward  HR 
implementation.  As  an  example  consider  a  simple  1st  order  IIR  filter 


# 


4 


art  =  ai  ♦  -f  Ujt  (2.16) 

As  in  the  FIR  case  we  introduce  a  delay  as  follows. 

Xk  =  D’^{ai*Xk~i+Uk)  (2.17) 

Now  we  use  one  delay  to  form  the  real  multiplier: 

Xk  =  D^~^{aiD*Xk.i  + Duk)  (2.18) 

=  D^~'^{ai*Xk-\  +  Duk)  (2.19) 

We  use  one  more  D  to  form  the  real  adder 

Xk- D^~'^{aiD*Xk^i+Uk)  (2.20) 

Finally  assume  D  =  z~^  and  Xk-i  -  z~^Xk  then 

Xk  -  z''^'^^{aiiz~'^Xk+z~^Uk)  (2.21) 

,=  z'‘^'^\aiixk+Uk)  (2.22) 


For  Xk  —  Xk^  N  must  be  zero.  However  iV  >  1  is  required  for  a  causal  filter. 

Solutions  to  this  problem  exist.  They  use  a  higher  order  difference  equation 
representation  of  the  filter  with  equivalent  characteristics,  so  that  feedback  loop  delay 
greater  than  1  can  be  tolerated.  More  details  can  be  found  in  the  literature  [Ref.  14] 
regarding  pipelining  of  IIR  digital  filters. 
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B.  SIGNAL  CHARACTERIZATION 


The  signal  yk  of  equation  (1.1)  can  be  characterized  either  in  the  time  domain 
or  frequency  domain  and  the  signal  czin  be  easily  transformed  from  one  domain  to  the 
other  via  a  linear  transformation. 

1.  Time  Domain  Representation 

In  the  time  domain  representation,  the  signal  yk  is  modelled  as  the  output 
of  a  linear  system  which  may  be  either  ein  all  pole  or  a  pole-zero  system  [Ref.  13]. 
The  input  is  assumed  to  be  white  noise  but  this  input  to  the  system  is  not  accessible 
and  is  only  conceptual  in  nature.  Let  the  system  under  consideration  be 

P  9 

yk  =  Yl  +  Z)  (2.23) 

i=l  j=0 

where  7*  is  a  white  noise  and  y*  is  the  output  of  the  system.  Equation  (2.23)  also 
can  be  represented  as  a  Transfer  Function(TF)  [Ref.  17] 


or 


where 


and 


(2.24) 

(2.25) 

A{z)  =  1 

»=i 

(2.26) 

B{z)  =  Y^bj+xz~^ 

(2.27) 

we  define  the  paiameter  vector  as 


p‘  =  [ai,a2,-**,ap,6i,62,**-,6, 


(2.28) 


The  parameter  p  completely  characterizes  the  signal  x*  The  above  system  defined 
by  the  equation  (2.23)  is  also  known  as  an  Auto* Regressive-Moving- Average(ARMA) 
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model  [Ref.  12], [Ref.  15].  It  is  conventional  to  denote  a  p  poles  and  q  zeros  system 
as  ARMA(p,q). 

2.  Frequency  Domain  representation 

The  signal  yk  of  equation  (2.23)  can  also  be  represented  as  a  vector 

Yjv  =  [yi*-*yiv]  (2.29) 

which  can  be  transformed  into  the  frequency  domain  by  the  DFT  relation 

=  (2-30) 

where  W  =  e~^  and  j  =  y/^.  In  general  y*  is  a  complex  quantity.  However  we 
restrict  our  interest  to  the  power  spectrunr  a  real  quantity  given  as 

»k=9kgl  (2.31) 

The  series  5^ 

Sn  =  [«i  •  •  •  «/v]  (2.32) 

is  another  form  of  representing  the  signal  and  is  designated  as  the  power  spectrum  of 
the  signal  z*.  Even  though  the  signal  in  this  domain  is  very  convenient  to  handle, 
computational  complexity  and  other  considerations  limit  its  use  for  online  application. 
Also  in  this  representation  phase  information  of  the  signal  is  lost. 

C.  SPECTRUM  ESTIMATION 
1.  Least  Square  Estimator 

Consider  the  signal  yjt  the  output  of  a  linear  system  excited  by  a  white 
noise  as  given  by  the  equation  (2.23).  Now  our  problem  is  to  estimate  the  parameter 
vector  p  by  obtaining  successive  measurements  of  y*.  We  shall  define  the  vector 

^k  ~  *  *  *  yVk—p)  *  ‘  *  >7^-9)  (2.33) 
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and  Zk  —  Vk  so  that  the  difference  equation  (2.23)  can  be  rewritten  as 


2k  =  Pk^k 

an  elegant  solution  [Ref.  2]  for  the  above  problem  can  be  written  as 

Pi  =  Pi-1  +  HiXtCi 


(2.34) 


(2.35) 


where 


(2.36) 


ek^zk-  p*.iXi 

The  block  solution  for  the  equation  (2.34)  can  be  also  written  in  the  form 


Pi  =  Hi  Xi^i^ 


(2.37) 


(2.38) 


Equation  (2.37)  can  also  viewed  as  Axk  -  B'^k  and  this  is  represented  in  Fig  2.3. 

Computing  the  matrix  Hi  via  equation  (2.36)  is  computationally  inefficient 
and  a  recursive  solution  via  the  matrix  inversion  lema  (Ref.  3]  is  preferred  and  is  given 


Hi  -  Hi-(  - 


Hi-jXixiHi-i 


(2.39) 


ll+xiHi-,XiJ 

Staticians  refer  to  the  Hi  matrix  as  the  eeuertonce  matrix  while  optiaii2ation  p  ople 
refer  to  Hi  as  the  //esston  of  the  objective  function,  where  the  objective  function  is 
given  by 

A  =  (2.40) 

Appendix  A  gives  a  computer  program  based  on  this  approach  that  was  used  for  this 
thesis  work  to  identify  system  parameters  via  Equations  (2.35),  (2.36),  (2.37),  (2.39). 
The  above  method  is  known  to  work  well  as  long  as  the  measured  value  of  the  slate 
of  the  filter  is  free  from  noise.  In  the  event  Zk  is  not  noise  free,  estimates  are  known 


to  be  biased. 
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2.  Maximum  Likelihood(ML)  Estimator 


The  ML  estimator  for  the  above  problem  is  conventionally  obtained  by 
considering  a  log-likelihood  function  (Ref.  16]  and  minimizing  it.  However  for  getting 
a  better  physical  understanding  the  approach  given  by  Tottn^lRef.  2]  will  be  followed. 
Consider  a  system  as 


Vk  = 


B(z) 


Aiz) 


'Ik 


and  an  auxiliary  system  referred  to  as  the  model  is  given  as 


(2.41) 


ik  = 


(2.42) 


and  we  dehne  the  function  to  be  minimized  as  L  =  clip)  where  Ck  yk  —  Xk. 
To  minimize  this  function  the  derivatives  and  the  Hessian  of  the  function  L  where 
i  —  (yfc  ~  Xk)^  need  to  be  obtained.  First  take  partial  derivative  of  L  with  respect  to 
each  paranreter.  For  example 


=  2(y*  -  «)l-gg;-)  = 


(2«) 


The  partial  ^  can  be  obtained  by  diflerentiating  equaUon(2.42)  as 


±(\ . +  ],V 

{2.44) 

U(i)  J 

dbi  “*  OiS“’  -  **  032“^  j  ) 

(2.45) 

Similarly  the  other  partial  ^ 

cat)  be  obtained  as 

d  \ 

■  1  MWI’ 

(2.46) 

(2.47) 
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since  =  0  as  B{z)  has  no  terms  containing  Oi.  Using  equations  (2.45)  and 

(2.47)  the  partials  axe: 

oai  L[>i(^)rJ 

It  = 

00,1  lA{z)] 

(2.51) 


by  defining 


Pfc  =  (<11502,03,61,^2] 

—  (^fc-l.^»--2i^’*-3iOjb,  Wfc-ll 

~  lJ/fc-Hyfc-2)  yk-3>  Ofc, 


the  Hessian  may  be  approxiniated  by: 


(Ex.'*|) 

\»sl  / 


Using  the  above  equations  the  parameter  can  be  estimated  via 

Pfc  =  Pfc-i  +  /iHfeX*efc 


(2.52) 

(2.53) 

(2.54) 


(2.55) 


(2.56) 


where  p  is  the  step  size  in  incrementing  the  parameter  vector.  A  block  schematic  for 
the  algorithm  was  given  in  Fig  2.4.  A  serious  drawback  of  this  method  is  the  stability 
while  incrementing  the  parameters. 

3.  Other  Methods 

There  are  many  other  methods  such  Cv.**  the  Instrumental  Variable  (IV) 
approach,  the  Approximate- Maximum- Likelihood  (AML)  Method  and  a  combination 
method  IV- AML  available  in  the  literature  (Ref.  2]. 
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D.  SELECTION  OF  ESTIMATION  TECHNIQUE 

Before  we  select  the  appropriate  method  for  the  problem  of  interest,  model-order 
problems  must  be  considered. 

1.  Model  Order  Problems 

In  the  Least-Square  Estimator  for  the  parameter  vector  p  of  the  signal 
Xk  we  need  to  assume  the  dimension  of  p  or  the  order  of  the  system.  Since  the 
system  order  is  often  not  known  a  higher-order  model  estimate  is  assumed.  Now  we 
investigate  higher-order  model  fits  for  a  lower-order  data. 

Consider  the  system  equation  (2.23)  with  parameter  given  as 

p‘  =--  (1.7,  -1.53, 0.648, 1, 0.6]  (2.57) 

For  the  above  system  a  Random  Binary  Noi8e(RBN)  was  given  as  the  input  u*  and 
the  corresponding  output  Xk  is  obtained.  For  the  time  series  u*  and  tk  an  ARMA(6,1) 
model  was  fit  using  the  program  given  in  the  appendix  and  the  parameter  was  esti¬ 
mated  as 

p‘ «  [2.26533,-2.79577,2.17249,-1.09877,0.446598,-0.110151,0.998175]  (2.58) 

Fig  2.5  gives  the  parametric  spectrum  for  the  above. 

For  the  same  data  an  ARMA(4,1)  model  was  used  and  generating  a  pa¬ 
rameter  vector 

p‘  =  [2.16988,-2.41286,1.47436,-0.365412,0.999317]  (2.59) 

Fig  2.6  gives  the  parametric  spectrum  for  the  above  values. 

Lastly  for  the  same  data  an  ARMA(8,1)  model  was  used  and  generating  a 
parameter  vector 

p‘  =  [2.276, -2.8'’5, 2  267,-1.282,0  703,-0.353,0.143,-0.0035,0.9999]  (2.60) 
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Fig  2.7  gives  the  corresponding  spectrum. 

It  is  interesting  to  compare  these  with  the  actual  parametric  spectrum  given 
in  Fig  2.8  For  the  sake  of  completion  and  comparison,  the  spectrum  given  in  Fig  2.9 
of  the  output  time  series  Xk  computed  via  a  DFT  program  is  also  included.  This 
simulation  demonstrates  that  the  existence  of  multiple  solutions,  hence  an  important 
step  in  the  estimation  procedure  is  to  choose  an  appropriate  model  order. 

2.  Choice  of  the  Method 

Since  the  filter  hardware  must  be  implemented  in  real-time  at  the  frequen¬ 
cies  of  10  to  20  Mhz  computational  complexities  must  be  kept  to  a  minimum.  These 
methods  although  providing  good  performance,  are  not  well  suited  to  hardware  im¬ 
plementations.  This  motivates  looking  into  new  filtering  methods. 


Pig  2.7  Fig 
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Fig 


III.  METHODS  for  FILTERING 


A.  FILTERING  TECHNIQUES 

Filtering  in  a  broad  sense  is  selectively  suppressing  a  portion  of  the  spectrum 
of  the  given  signal.  In  this  chapter  several  filtering  techniques  are  explored  in  the 
light  of  the  narrow-band  interference  problem  in  order  to  identify  applicable  filtering 
techniques. 

1.  Filtering  via  FFT 

Any  given  signal  can  be  conveniently  transformed  into  its  power  spectrum 
via  equation  2.31.  Assuming  the  spectrum  of  the  desired  filter  to  be  defined  by  Wk- 
The  weighted  output  is  given  by: 

Vk  =  SkWk  (3.1) 

and  the  corresponding  filtered  output  yu  is  obtained  by  taking  the  inverse  DFT  of  the 
signal  Vk>  A  simple  block  schematic  representing  this  idea  is  given  in  Fig  3.1.  Details 
of  this  approach  are  available  in  various  references  [Ref.  22]. 

Fast  algorithms  such  as  the  FFT  for  computing  the  DFT  make  it  possible 
to  do  the  above  process  in  real  time  by  using  dedicated  hardware.  Honeywell  makes 
the  HDSP66110  and  HDSP66210  Digital  Signal  Processing  chip  pair  which  are  ideally 
suited  for  these  applications.  This  chip  pair  can  perform  a  single  complex  multiply  in 
40n8  [Ref.  28).  However  the  problem  of  filtering  adaptively  [Ref.  10]  (i.e  varying  Wk 
according  to  a  criterion)  demands  more  computing  power  which  can  be  obtained  by 
augmenting  the  DSP  chip  pair  with  a  processor  like  the  R3000  which  has  a  Reduced 
Instruction  Set  Computer(RISC)  architecture  with  a  high  instruction  execution  rate 
of  approximately  25  Million  Instructions  Per  Second(MIPS)  [Ref.  27]. 


23 


2.  Recursive  DFT 

Implementation  of  a  higher  order  FIR  filter  using  a  recursive  DFT  [Ref. 
23], [Ref.  24]  is  also  very  convenient  for  eliminating  narrowband  interference.  The 
basic  concept  is  that  the  FIR  filter  is  expressed  as  a  product  of  two  filter  sections. 
One  section  is  a  filter  with  its  zeros  being  equally  spaced  on  the  unit  circle.  This 
is  achieved  by  a  delay.  The  second  section  is  a  pole-producing  section.  Pole-zero 
cancellation  results  in  the  desired  FIR  filter.  More  details  can  be  seen  in  the  references 
[Ref.  23], [Ref.  24]. 

3.  Adaptive  Filtering 

Adaptive  filters  can  be  placed  into  four  classes  based  upon  the  choice  of 
the  training  sequence  and  the  reference  model  for  adaptation.  Simon  [Ref.  21]has 
classified  the  Adaptive  Filters  into  the  following  four  classes 

•  Identification  -  Class  I 

•  Inverse  Modelling  -  Class  II 

•  Prediction  -  Class  III 

•  Interference  Canceling  -  Class  IV 

Fig  [3.2  -  3.5]  give  the  block  schematics  of  the  various  classes  of  adaptive 
filters.  In  our  problem,  the  reference  model  is  naturally  a  narrow-band  bandpass  filter 
since  the  interference  signal  is  a  narrow-band  signal.  The  class  of  filtering  that  best 
suits  our  problem  is  a  combination  of  Class  III  and  Class  IV  type  of  filter  [Fig  3.4 
and  3.5].  This  logically  points  to  an  adaptive  notch  filter.  Due  to  apriori  knowledge 
of  the  interference  signal,  only  the  parametric  approach  was.  However  the  DFT- 
based  techniques  may  also  offer  a  good  solution  and  should  be  the  subject  of  future 
investigations. 
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Identification  (  Class 


Delay 


B.  ADAPTIVE  NOTCH  FILTERS 


Notch  filters  for  removing  multiple  narrow-band  interference  can  be  categorized 
into  four  broad  categories  illustrated  in  Fig  3.6  through  3.9.  The  first  two  categories, 
Figures  3.6  and  3.7,  are  cascaded  second  order  notches  with  each  second-order  section 
removing  one  frequency.  The  next  two  categories, 

Figures  3.8  and  3.9,  are  higher-order  notches  that  eliminate  multiple  frequen¬ 
cies.  In  all  of  the  categories,  it  is  possible  to  use  FIR  filters  (t.e  all  zero  filters)  which 
are  easily  pipelined  and  can  be  made  truly  linear  phase.  However,  HR  notch  filters 
require  substantially  fewer  multipliers  and  adders  than  FIR  notch  filters.  Thus  HR 
pipelining  may  become  an  important  issue.  In  this  thesis  we  limit  our  discussion  only 
to  the  first  two  categories  illustrated  in  Figures  3.6  and  3.7. 

1.  Second-Order  Cascaded  Notch  Filters 

The  second-order  notch  filter  is  used  in  cascade  and  in-line  with  the  signal 
as  shown  in  Figure  3.6.  The  transfer  function  for  such  a  notch  filter  is  given  by  Kwan 
and  Martin  [Ref.  8]  as: 


1  -  m,{z) 


ki  (l-fz-‘)(l-z->) 

2  1  -  (2  -  fca  -  -f  (1  -  iba)2-a 


(3.2) 

(3.3) 


_  2  -  l-i  1  -  +  5-» 

“  2  1  -  (2  -  t, -tf)j-‘ +(!-*, )j-> 

For  arbitrary  values  of  ki  and  kt ,  this  is  a  symmetric  notch  filter  with  unity 
gain  at  DC  and  the  Nyquist  frequency.  If  itj  is  kept  constant,  then  the  3db  notch 
width  is  also  kept  constant.  Thus  ki  may  be  adapted  to  remove  one  narrow-band 
signal.  A  cascade  of  such  filters  can  be  used  to  remove  multiple  narrow-  band  signals 
[Ref.  8].  Constants  I'j  and  1*3  are  related  to  the  pole  radius  r  and  the  normalized  pole 
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Cascade  of  Second-Order  Notch  Filters 
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Higher-Order  Notch  Filter 


N-th  order 
M-Bandpass 
Filter 


frequency  Bp  as 


$ 


A 


=  (1  —  r^)  (3.5) 

=  \/l  +  >  ^  -  2rcosBp  (3.6) 

It  is  important  to  bear  ic  mind  that  the  above  transfer  function  has  unity  g^Lin  at 


Bpeak  =  2arcsin 


(3.7) 


which  is  different  from  Bp  [Ref.  8] 

2.  Second-Order  Cascaded  Signal  Canceler 

The  cascaded  second-order  signal  canceler  approach  shown  in  Figure  3.7 
has  the  advantage  that  the  desired  signal  does  not  pass  through  the  adaptive  filter, 
instead,  the  band-pass  filter  is  used  to  detect  the  narrow-band  signal  which  is  then 
subtracted  from  the  desired  signal.  A  constant  3db  bandwidth  notch  can  be  achieved 
by  selecting  a  band-pass  filte"  with  the  tramsfer  function: 


Htpiz) 


^2(1  - 
2D(z) 

D{z) 


(3.8) 

(3.9) 


where  D{z}  is  the  same  denonrinator  as  Hn{z)  in  equation  (3.4)  The  signal-canceler 
structure  is  also  nice  for  adaptation  because  it  is  relatively  easy  to  generate  sensitivity 
functions  which  are  related  to  the  gradient  of  Hbp{z)  with  respect  to  the  frequency 
parameter  fcj.  Figure  3.10  shows  the  block  diagram  of  an  adaptive  version  of  this 
filter.  The  sensitivity  function  s{n)  is  obtained  by  differentiating  the  error  signal 
c(n)  with  respect  to  the  parameter  fcj.  This  can  be  easily  derived  by  noting  that 
e(n)  =  x(n)  -  j/(n)  and 


de{n)  _  dy{n) 
dki  dki 


(3.10) 
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to  get  the  above  partial  derivative  we  consider 

y{n)  =  Hbp{z)x{n) 

m. 


D{z) 


i(n) 


(3.11) 

(3.12) 


the  above  equation  (3.12)  can  be  easily  differentiated  with  respect  to  ki  by  using  the 
equation  (2.47)  to  obtain: 


dy{n) 

dki 


■^2hz-' 


;(n) 


(3.13) 


by  recognizing  that  Hbp{z)  =  from  equation  (3.12)  the  above  equation  can  be 


written  as 


and  by  defining 


dyjn) 

dki 


Hb,{z) 


2kiz' 


H»iz)  = 


Diz) 
2kiz~^ 

~WT 


a:(n) 


(3.14) 


(3.15) 


we  get  the  sensitivity  function  as 


>(n)  =  ^  =  H.(z)  = 


(3.16) 


The  equation  (3.16)  is  given  &s  a  block  schematic  in  Fig  3.10.  The  parameter  ki  may 
then  be  adapted  by  the  formula: 


ki{n  +  1)  =  fci(n)  -  /ie(n)- 


s(n) 


a.  Computing  |l5(n)|l* 

Ideally  ll5(n)|l^  can  be  calculated  by  the  relation 


(3.17) 


lk(«)ll*=  E  ^(0^ 

is:n-N 


(3.18) 
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The  above  scheme  computes  the  average  of  the  past  N  samples.  However  a  weighted 
average  of  the  past  samples  with  highest  weight  for  the  recent  sample  can  be  done 
quite  easily  by  a  first-order  low  pass  filtering  [s(n)]''^  as  follows: 

Vn  =  Vn-iA  -I-  {1  -  A)s^{n)  (3.19) 


However  the  above  lowpeiss  filtering  can  also  be  carried  out  by  a  second  order  filter 
such  as: 


6i(1-22-'  -f  z-2) 


Vn  = 


[  1  -  2rfkz-^  -I-  J 

6i(l  —  2cos{h6)z~^  -f  z~'^)^ 


1  —  2r jkz~'^  -f  z~^ 
Yet  another  way  to  estimate  lls(n)|l^  i  simply: 


[s(n)]^ 


N 


Vn  =  S  Win)] 


(3.20) 

(3.21) 


(3.22) 


f=i 

the  above  form  is  defined  as  zero  order  forgetting.  It  should  be  noted  that  equation 
(3.21)  was  used  by  Kwan  and  Martin  (Ref.  8]. 


C.  MULTIPLE  NARROWBAND  SIGNAL  SUPPRESSION 

The  second-order  implementations  of  section  (1.)  and  (2.)  offer  considerable 
advantage  both  in  hardware  complexity  and  in  adaptive  performance  ,for  multiple 
narrow  band  signal  suppression  [Ref.  8], [Ref-  25], [Ref.  26). 

1.  The  Kwan  and  Martin  Filter 

In  a  recent  paper  by  Kwan  and  Martin  [Ref.  8],  the  problem  of  detecting 
and  enhancing  sinusoidal  signals  in  the  presence  of  noise  is  addressed  with  a  cascade 
of  HR  adaptive  notch  filters  which  are  used  to  eliminate  the  sinusoids.  Each  of  the 
sinusoids  is  eliminated  by  a  bandpass  filter  whose  output  is  an  enhanced  version  of 
one  of  the  sinusoids.  Hence  this  remarkable  structure  can  perform  both  tasks  with  a 
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single  adaptive  filter  configuration  which  is  shown  to  be  highly  robust  and  performs 
extremely  well. 

The  major  disadvantage  of  the  Kwan  and  Martin  structure  is  that  the 
number  of  biquad  sections  needed  in  the  adaptive  filter  configuration  is  given  by 
N(N-f3)/2,  where  N  is  the  number  of  sinusoids  to  be  detected  and  removed.  This 
becomes  impractical  in  real-time  situations  with  more  than  4  sinusoids  due  to  the 
geometric  increase  in  the  required  hardware. 

a.  Kwan  and  Martin  Structure 

The  Kwan  and  Martin  structure  consists  of  a  cascade  of  HR  notch 
filters  one  stage  of  which  is  shown  in  Figure  3.11.  Each  stage  consists  of  a  bandpass 
filter  with  zeros  at  DC  and  the  Nyquist  frequency  and  unity  gain  at  its  peak  frequency 
a>i.  Such  a  filter  would  have  the  following  z-domadn  transfer  function 


where 

n  =  pole  radius  of  the  i-th  section 

oji  =  peak  frequency  of  the  i-th  section 

Ut  =  sampling  frequency 

Kwan  and  Martin  identify  two  different  methods  for  adapting  the 
filter.  Most  of  their  derivation  is  based  on  what  they  call  the  constant  bandwidth  filter 
in  which  the  pole  radius  r<  is  a  constant  and  only  the  frequency  w,'  is  adapted.  An 
alternative  approach  which  keeps  a  constant  Q  is  also  discussed  in  Kwan  and  Martin 
paper.  In  addition,  Kwan  and  Martin  select  the  adaptive  quantity  in  such  a  way  that 
it  is  fairly  easy  to  determine  the  notch  frequency  from  the  adaptive  parameter.  FVom 
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Figure  3.10,  we  see  that  the  notch  filter  for  each  section  is  the  difference  between  1 
and  the  bandpass  filter,  hence: 


=  1  -  K(^)  {3-24) 

b.  Calculation  of  the  gradient 

The  basic  structure  of  the  Kwan  and  Martin  adaptive  filter  shown  in 
Figure  3.12  is  a  cascade  of  N  sections  of  the  form  of  Figure  3.7,  The  overall  transfer 
function  is  given  by 


T{z) 


N 


n  (i  -  K) 

z _ « 


(3.25) 

(3.26) 


Kwan  and  Martin  choose  as  their  objective  function  J (2)  the  square 
of  the  output  of  the  final  stage  of  the  cascade: 


Jiz)  =  {E,{z)Y 

=  ir(^)f[A'(.)]» 


(3.27) 

(3.28) 


Hence,  the  gradient  of  the  objective  function  J{z)  is  given  by 


(3.29) 


Thus  in  order  to  find  the  gradieht  with  respect  to  the  adaptive  pa¬ 
rameters  kj  ,  we  must  take  the  partial  derivative  of  T{z)  with  respect  to  each  kj 


dT{z)  _  f  fMjz) 


(3.30) 
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From  equation  (3.24)  we  have 


^  a[l  -  Hi,{z)] 

dkj  dkj 

^  dHUz) 

dkj 


(3.31) 

(3.32) 


Using  the  rule  for  differentiation  as  given  by  the  equations  (2.45)  and  (2.47)  and  from 
equation  (3.23)  we  have 


dkj  “  ^  D{z) 

=  HUz)HUz) 


(3.33) 

(3.34) 


where  D{z)  is  the  denominator  of  Hi,  is  given  by  the  equation  (3.15). 

Substituting  equation  (3.34)  into  equation  (3.30)  we  obtain  the  over 
all  sensitivity  for  the  jih.  ki  parameter  as 


s 


iml 


i-i 


n ««(») 


s 


By  recognizing  that 

=  I  n 

and  by  substituting  equation  (3.37)  in  the  equation  (3.36)  we  get 


(3.35) 

(3.36) 

(3.37) 


mz) 

dkj  ~ 


ji“3 

ial 


«Liz) 


(3.38) 


Figure  3.12  shows  the  Kwan  and  Martin  realization  of  the  complete 
adaptive  system.  The  difficulty  in  the  Kwan  and  Martin  approach  is  in  the  generating 
of  the  product  of  notch  filters  without  the  notch  filter  j,  as  required  in  equation  (3.36). 
To  generate  this  product  for  each  section,  would  require  -  I  biquads  per  section 
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resulting  in  a  total  of  N{N  —  1)  biquads  just  to  generate  the  product.  Kwan  and 
Martin  are  able  to  reduce  this  by  uring  the  output  of  the  bandpass  filter  as  the  input 
to  their  cascade  via  equation  (3.37).  Since  this  output  already  has  (j  —  1)  of  the 
required  Hf^r{z)  factors  in  it,  only  (N—j)  additional  biquads  are  needed  for  a  total  of 
0.5(iV^  -  N).  Adding  this  to  the  N  biquads  required  to  realize  the  cascade  of  notch 
filters  and  the  N  biquads  required  to  realize  the  Hl{z)  factors,  yields  a  total  number 
of  biquads  given  as  0.5(A'^  +  ‘i)N. 

2.  The  New  Structure 

Figure  3.13  shows  the  improved  adaptive  notch  filter  structure  [Ref.  25), [Ref. 
26).  The  key  to  the  improvement  is  the  recognition  that  the  output  E\{z)  —  T{z)X{z) 
for  the  cascade  of  the  notch  filters  can  be  written  both  as  a  product  of  the  individual 
notch  filter  section  transfer  functions  times  the  input  and  in  terms  of  the  input 
X{z)  minus  the  outputs  Y'is)  of  the  bandpass  filters 


T(r)A'(,) 

(3.39) 

(3.40) 

Jia)  J 

AW  -  (zno) 

(3.41) 

To  gel  the  product  of  without  the  term  i  =:  ji  as  required  in  the 
equation  (3.36),  we  may  use  equation  (3.41)  to  simply  add  back  the  term  Y^is) 


A'(r) 

-  A'(i)- 

u;;  j 

«  /i:t(2)+y^(5) 


(3.42) 

(3.43) 


Figure  3.13  makes  use  of  this  fact  to  generate  the  gradient  needed  for  the 
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New  Structure 


Fig  3.13 


adaptive  process.  From  the  Figure  3.13,  we  can  see  that  the  total  number  of  biquads 
required  is  N  for  the  cascade  of  notch  filters  plus  2N  for  the  Hlp{z)Hl{z)  required  for 
adaptation,  minus  1  at  the  last  stage,  since  the  last  stage  does  not  need  the  extra 
Hlj,{z)).  Thus  we  have  the  number  of  biquads  required  in  the  new  structure  is  (3N-1). 


* 


* 
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IV.  MODELLING  and  SIMULATION 


A.  MODELLING  OF  SIGNALS 

In  order  to  test  various  algorithms  and  to  evaluate  their  probable  performance 
in  the  real  environment,  it  was  necessary  to  develop  a  meaningful  simulation  of  the 
real  situation.  To  achieve  this,  the  following  four  testing  categories  were  developed: 

(i)  Sinusoidal  signals  with  white  gaussian  noise 

(ii)  Narrow  Band  Noise  with  white  gaussian  noise 

(iii)  Bi-Phase  Shift  Keying  (BPSK)  sequence 

(iv)  Frequency  Shift  Keying  (FSK)  sequence 

1.  Sinusoidal  Signals 

In  order  to  generate  sinusoidal  signals  with  minimum  computational  burden 
placed  at  different  frequencies  a  second  order  AR  process 

X*  ss  -  4-8  (4.1) 

was  used.  Initial  conditions  are  very  important  and  they  arc  chosen  such  that  r_i  =  0 
and  s  giving  a  unit  amplitude  sinusoidal  signal.  The  9i  value  is  between 

0  to  160  and  n  is  the  number  of  frequencies  desired.  The  required  signal  needed 
to  input  into  the  adaptive  algorithm  is  given  as 

fl 

y*  “514  + 7k  (4.2) 

ia>l 

where  %  is  a  white  gaussian  noise  Af(0,£r’). 

2.  Narrow  Band  Noise 

The  narrow  band  noise  signal  is  generated  using  the  difference  equation 

4  =  2rcos(«i)xi.,  -  r’lLa  +  4  (4.3) 
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where  6i  decides  the  placement  of  the  noise  in  the  spectrum  and  r  controls  the  band¬ 
width  of  the  noise.  The  u\  is  simply  a  uniformly  distributed  noise  taken  at  different 
instants.  The  desired  signal  yt  is  obt^ned  via 

+  (4.4) 

i=l 

3.  Bi-Phase  Shift  Keying  Sequence 

The  generation  of  BPSK  signal  has  three  distinct  three  parts; 

(i)  Generation  of  Random  Binary  Sequence(RBS)  is  achieved  by  passing 
uniformly  distributed  noise  through  a  hardlimiter  (An  important  note  is  that  the 
interval  between  the  two  consecutive  bits  of  RBS  is  j-); 

(ii)  Another  sequence  biniary  numbers  is  a  spreading  code  or  sequence. 
The  specific  sequence  used  in  a  given  communications  system  is  normally  not  available 
to  anyone  but  to  the  designated  receiver,  (In  this  particular  simulation  we  have 
generated  the  spreading  aequeuce  by  passing  uniformly  distributed  noise  through  a 
hardlimiter.  The  bit  interval  is  ^); 

(iii)  Phase  eacodiug  (i.e.)  mapping  the  given  binary  signal  which  is  the 
exclusive  or  of  (i)  and  (ii)  as  0  orir  at  appropriate  santple  time. 

The  output  of  the  first  hardlimiter  is  storfid  in  an  array  z.  Output  of  the 
second  hardlimiter  is  stored  in  the  array  y  This  information  is  relric\^  by  a  subtle 
use  of  the  array  index  given  as  i  e  Jb/»  where  i  is  the  index,  k  is  the  discrete  sample 
number  and  is  the  bit  rale  of  the  intelligence.  Similarly  another  index  ;  is  generated 
using  j  fi=  kfe  where  /«  is  the  chip  frequency.  Generation  of  the  indices  is  the  key 
thing  in  this  simulation.  The  desired  signal  now  is  given  by  the  equations 


tUt  =  Acos{2irfok  -f 

(4.5) 

(4.6) 
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P=1 

where  i  and  j  are  indices  of  the  arrays  as  defined  earlier.  It  is  an  important  point 
to  note  that  the  characteristic  of  the  signal  wl  are  defined  by  the  parameter  p  — 

{/Oi/o/b}- 

This  signal  is  not  really  a  simple  BPSK  signal  but  it  has  an  additional  fea¬ 
ture  of  spreading  the  spectrum  by  controlling  the  chip-frequency  and  carrier  frequency 
and  information  rate.  A  block  diagram  of  the  scheme  is  given  in  Fig  4.1. 

The  desired  signal  y*  is  given  by; 

y*  =  oti  +  0k  (4.8) 

where  ak  is  the  set  of  narrowband  BPSK  signals  placed  at  different  places  in  the 
frequency  spectrum  and  0k  is  the  broad  band  BPSK  signal  generated  for  a  specific  p 
value, 

4.  CVequency  Shift  Keying  Sequence 

Generating  this  sequence  needs  a  random  binary  intelligence  signal.  This 
was  once  again  is  achieved  by  passing  a  uniformly  distributed  noise  through  a  hardlinv 
iter.  The  output  of  this  hard  limiter  stored  in  an  array  t.  An  index  t  is  chosen  such 
that  t  «  kfk  where  /s  is  the  haud-rate  of  the  information  and  k  is  discrete  sample 
number.  Now  the  desired  signal  is  generated  via 

Sk  ®  2cos(0k)sk^i  ~  H-7  (4.0) 

=  (4.10) 

ya  =  ak  +  7i  (4.11) 

where  $  is  the  carrier  frequency  and  6  is  the  depth  of  the  frequency  modulation. 
Initial  conditions  are  very  important  and  they  are  chosen  such  that  s-i  =  0  and 
s-3  =3  -Mn{6)  giving  an  unit  amplitude  sinusoidal  signal. 
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B.  SIMULATION 


The  adaptive  digital  filter  algorithm  as  described  in  the  earlier  chapters  is  sim¬ 
ulated  and  tested  using  synthetic  data.  Simulation  was  carried  out  on  a  VAX  11/785 
computer  using  Fortran  77.  Listing  of  the  program  used  is  enclosed  in  the  appendix. 
The  adaptive  filter  parameters  of  interest  are 

(a)  Sharpness  of  the  notch  filter  defined  by  pole  position  (r^  =  1  —  K2) 

(b)  Step  size  in  the  parameter  update  procedure  (/r) 

(c)  Time  constant  of  the  fading  filter  (A)  (  refer  to  equation  3.19) 

(d)  Model  order  (A'=  #  of  2nd  order  filters) 

(e)  Order  of  the  incoming  signal  or  number  of  interferers  (m) 


1.  Kwan  &  Martin  Algorithm 

In  simulating  this  algorithm,  the  most  important  thing  is  the  implemen¬ 
tation  of  the  notch  filter.  Fig  4.2  gives  the  structure  of  the  algorithm.  Let  x{n)  be 
the  input  to  the  Adaptive  Filter  and  e(n)  the  desired  output.  This  desired  output  is 
obtained  by  passing  the  input  x(n)  via  a  cascade  of  N  notch  filters  as 

=  (4.12) 

Output  at  the  intermediate  jth  section  of  the  band-pass  filter  is  designated  as  y^{n) 
as  shown  in  Fig  4.2  Then  the  sensitivity  of  the  k\  parameter,  is  obtained  by  passing 
this  j/^{n)  through  another  cascade  of  (j  -  1)  notch  filters  given  as 

An)  =  ■  ■  ■  ilh«L}An)  (4.13) 

The  filter  transfer  functions  Hj^{z)  and  Hi,  were  defined  in  the  earlier  chapter  by 
equations  (2.21)  and  (2.15). 
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2.  The  New  Algorithm 

The  primary  diflFerence  between  the  New  Algorithm  and  the  Kwan  and 
Martin  algorithm  is  in  the  method  of  obtaining  the  sensitivity  of  the  jth  coefficient. 
This  difference  can  be  seen  in  the  Fig  4.3.  The  output  y^{n)  at  the  jth  notch  Fig  4.3. 
filter  is  added  to  the  over-all  output  e(n)  and  this  summed  output  is  passed  through 
a  cascade  of  only  two  filters  as  shown  in  Fig  4.3  and  can  be  given  as 

«>(n)  =  {e(n)  +  !,>(n)}  (4.14) 

3.  Forgetting  Filters 

For  parameter  incrementation  via  equation  (3.17)  we  need  to  obtain  ||s'^(n)l|. 
This  could  be  achieved  by 

a)  Zero  order  forgetting  using  equation  (3.22) 

b)  1st  order  forgetting  using  equation  (3.19) 

c)  2nd  order  forgetting  using  equation  (3.21) 

4.  Stability 

The  parameter  incrementation  given  by  the  equation  (3.17)  can  be  recast 
by  using  either  of  the  forgetting  filters  given  above  as 

k\  =  k\  -  /i[v']~^s’(n)e(n)  (4*15) 

It  is  very  important  to  note  that  the  above  equation  (4.15)  has  a  close  resemblance 
with  the  ML  Estimator.  In  the  MLE  case  the  stability  while  incrementing  the  param¬ 
eter  is  an  important  factor.  A  similar  problem  exists  in  the  current  incrementation 
procedure  but  the  problem  is  solved  by  a  suitable  choice  of  parameters  and  modify¬ 
ing  the  incrementation  procedure  by  adding  an  additional  factor  to  equation  (4.15) 
as  follows: 

(4.16) 


k\  =  k\  -  ^e{n)s*{n) 


1  +  -^ 
_ «mo« 


+  V‘ 
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Fig  4.3 


This  modification  also  protects  against  possible  overflow  or  underflow  while  computing 
[u*]~\  In  addition  simulations  indicate  that  in  the  case  of  zero-order  forgetting,  this 
factor  is  required  in  order  to  obtain  convergence.  When  v^in  and  Vmax  are  zero  and 
infinity  respectively,  equation  (4.16)  reduces  to  equation  (4.15). 

In  this  incrementation  procedure  two  forms  of  the  denominator  polynomial 
of  the  Hbp{z)  were  considered: 


and 


D\z)  =  1  -  (2  -  1:2  -  hl)z~^  +  (1  “ 

(4.17) 

£>‘(2)  =  l-2K‘r2-^  +  rV2 

(4.18) 

where  k\  =  cos(0'). 

While  using  the  incrementation  procedure,  under  transient  conditions  poles 
of  the  D*{z)  cross  the  unit  circle  causing  instability  problems.  This  condition  was 
averted  by  checking  the  pole  position  after  incrementation  and  if  unstable  then  in¬ 
crementation  is  modified.  Checking  this  condition  for  D*{z)  given  by  equation  (4.17) 
calls  for  solving  a  quadratic  equation  at  every  parameter  increment  and  examining 
the  pole  position.  However  this  check  is  very  simple  for  D*{z)  given  in  equation 
(4.18),  since  we  need  only  to  check  k{  by  maintaining  l/cd  <  1.  Furthermore  the 
value  of  r  must  be  positive  and  less  than  unity  for  stability.  The  choice  of  the  r  is 
very  important  for  proper  fast  convergence.  Fig  4.4.  shows  the  frequency  response  of 
band-pass  filter  for  different  values  of  r  The  3dB  band- width  of  the  band  pass  filter 
is  related  to  r  and  is  given  [Ref.  29]  by 


Band  —  width  =  cos"^ 


(4.19) 


Note  that  under  limiting  conditions,  the  band-width  is  essentially  zero.  In  this  simu¬ 
lation  it  was  assumed  a  value  of  r  as  0.95. 
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A  large  number  of  simulations  [Ref.  29]  were  carried  out  using  sinusoidal 
inputs  with  and  without  noise  for  different  m  and  N  values.  In  these  simulations 
the  denominator  polynomial  used  was  given  by  equation  (4.18).  The  results  were 
tabulated  in  Table  3.1.  From  Table  3.1  it  is  seen  that  the  choice  of  the  step  size  is 
an  important  factor  and  step-size  must  be  optimized  for  a  specific  number  of  sections 
A.  The  values  of  v^in  and  Vmax  did  not  pose  any  problem  while  incrementing  using 
1st  order  or  2nd  order  forgetting  filters.  It  seems  that  a  1st  order  forgetting  filter  is 
more  effective  than  either  zero-order  or  a  2nd  order  filter. 

After  fixing  the  values  of  fi  and  A  the  algorithm  was  tested  using  Narrow 
Band  Noise  signal  y*  for  its  performance.  This  signal  was  used  only  to  tune  the  value 
of  r  (i.c.  sharpness  of  the  notch  filter).  The  following  simulation  results  axe  based  on 
parameter  adaptation  for  the  D*{z)  given  by  equation  (4.17). 

5.  Response  for  Sinusoidal  signal 

A  sample  simulatior  output  due  to  sinusoidal  signal  was  shown  in  Fig  [4.5 
•  4.7).  The  input  signal  is  composed  of  3  sinusoids  with  normalized  frequencies 

and  1^/,  and  WGN  with  =  1.  The  spectrum  of  this  signal  is  shown  in  Fig 
4.5.  After  passing  this  signal  through  the  adaptive  filter  we  could  see  that  interferers 
were  removed  and  only  the  noise  was  left  behind.  This  is  clearly  shown  in  Fig  4.6. 
The  adaptation  process  is  shown  in  Fig  4.7. 

6.  Response  for  Narrow  Band  Noise 

In  this  testing  category  sinusoids  were  replaced  by  narrow-band  signals. 
These  signals  were  generated  via  equation  (4.4).  Superimposed  on  this  signal,  WGN 
with  a  =  1  was  added.  A  typical  spectrum  of  this  signal  is  shown  in  Fig  4.8.  This 
signal  was  passed  through  the  adaptive  filter  and  these  narrow  band  interferers  are 
notched  out  and  only  WGN  is  left  out  as  shown  in  the  Fig  4.9.  Fig  4.10  shows  the 
corresponding  parameter  convergence. 
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7.  Tracking  FSK  Signal 

Transient  behavior  of  the  adaptive  filter  was  studied  by  applying  an  FSK 
signal  (generated  via  equation  4.11).  In  the  case  of  an  FSK  signal,  signal  spectrum 
constantlj  changes.  This  property  is  useful  for  testing  the  tracking  ability  of  the 
adaptive  filter.  In  fact  a  WGN  with  <t  =  1  was  also  added  to  the  signal.  An  adaptive 
filter  then  was  used  to  demodulate  the  signal,  by  tracking  the  spectrum.  This  is 
shown  in  Fig  4.13. 

8.  Suppressing  BPSK  Interference 

•  Having  seen  the  performance  of  the  adaptive  filter  under  various  conditions, 

it  is  only  needed  to  test  under  an  additional  constraint  of  broad  bsmd  noise.  In  this 
case,  the  broad-band  signal  was  swamped  by  narrow-band  interferers  and  WGN.  A 
typical  broad-band  signal  is  generated  by  using  a  BPSK  signal  generator  via  equation 
(4.8).  In  this  we  have  used  carrier  at  0.25  cycle/sample  and  the  chip  frequency  at  0.125 
cycles/sample.  The  narrow-band  interference  signals  are  also  generated  by  using  the 
same  BPSK  signal  generator,  but  with  different  parameters.  Interferers  are  chosen 
such  that  the  chip  frequency  is  fixed  at  cycles/sample,  while  carriers  were  chosen 
^/«»  l^/t  cycles/sample.  To  these  interferers  and  broad-band  signal, 

a  WGN  with  <r  =  1  was  added.  The  composite  signal  spectrum  is  shown  in  Fig  4.14, 
indicating  the  three  interferers,  broad-band  signal,  and  WGN.  This  signal  was  passed 
through  the  adaptive  filter.  Output  of  the  filter  is  shown  in  Fig  4.15.  In  this  figure 
it  is  clear  that  the  interferers  were  removed  by  the  filter  and  only  the  desired  signal 
►  was  left  behind. 

9.  Model  order  mismatch 

The  model  order  of  the  algorithm  was  fixed  at  3  while  BPSK  signals  were 
generated  with  2  interferers  t.e.  m  —  2  and  N  =  3.  Fig  3.16  shows  the  parameter  con- 
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vergeoce  under  these  conditions  and  Fig  3.19  shows  for  m  =  1  and  iV  =  3  conditions. 
We  could  notice  that  free  extra  paxameter(s)  wandering  due  to  the  mismatch.  This 
has  the  effect  of  notching  out  the  desired  signal.  This  can  be  solved  by  deliberately 
introducing  a  known  BPSK  signal  at  fixed  place. 
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Summary  of  Tests  on  the  New  Algorithm 


TABLE  4.1 


Case 

4  of 
sections 

Sines 

"  1 

^min 

Vnxax 

Noise 

zero  order 

1st  order 

2nd  order 

HSU 

iter 

mm 

iter 

mm 

iter 

la 

1 

1 

0.05 

20 

no 

0.05 

67 

0.05 

20 

0.05 

20 

Ian 

1 

1 

0.05 

20 

yes 

0.05 

WESSl 

0.05 

20 

0.05 

144 

lb 

1 

3 

0.05 

20 

no 

0.05 

0.05 

30 

25 

Ibn 

1 

3 

0.05 

20 

yes 

0.05 

jump 

0.05 

50 

"005^ 

100 

2a 

5 

3 

0.10 

Msm 

no 

0.017 

340 

T& 

200 

WT 

400 

2an 

5 

3 

0.10 

10 

yes 

0.017 

1200 

0.063 

400 

0.07 

450 

2b 

5 

5 

0.10 

10 

no 

0.017 

1410 

0.063 

150 

0.07 

500 

2bn 

5 

5 

0.10 

yes 

0.017 

3500 

175 

250 

2c 

5 

BIH 

0.10 

10 

no 

0.017 

1667 

0.063 

450 

IflGM 

BBSS 

2cn 

OH 

0.10 

10 

yes 

0.017 

700 

0.063 

167 

0.07 

KTff.ia 

3a 

7 

0.02 

50 

no 

0.02 

1667 

0.075 

500 

0.07 

1500 

3an 

10 

7 

0.02 

50 

yes 

0.02 

5000 

0.075 

750 

0.07 

1000 

3b 

10 

10 

0.02 

50 

no 

0.02 

7500 

0.075 

2000 

0.07 

5250 

3bn 

10 

10 

0.02 

50 

yes 

0.02 

1000 

0.075 

2000 

0.07 

2500 

•  In  all  cases  with  noise,  parameters  converge  to  exact  value  when  number  of 


notches  is  greater  or  equal  to  number  of  sinusoids.  Iterations  listed  in  table 
represent  convergence  to  three  decimal  places. 

•  In  all  cases  with  noise,  parameter  oscillate  around  the  exact  value.  The  number 
of  iterations  listed  in  the  table  is  the  number  of  iterations  required  to  establish 
this  oscillatory  pattern. 

•  In  some  cases  with  more  interfering  sine  waves  than  notches  will  jump  at  random 
between  sine  waves.  This  pattern  is  indicated  in  the  table  by  the  word  jump. 

•  Values  ol  and  Vmor  in  the  table  refer  only  to  zero-order  forgetting  which 
generally  will  not  converge  without  limits  on  w.  All  other  types' of  forgetting 
were  run  without  limits  on  w 
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V.  HARDWARE  IMPLEMENTATION 


Although  hardware  design  is  beyond  the  scope  of  this  thesis,  in  this  chapter  we 
shall  briefly  look  at  possible  hardware  configurations  with  a  view  toward  feasibility 
of  the  new  algorithm. 

A.  HARDWARE  FEASIBILITY 

In  recent  years  many  dedicated  digital  signal  processing  chips  have  become 
commercially  available.  Table  5.1  gives  characteristics  of  some  typical  chips  that  are 
currently  available. 


TABLE  5.1 


Feature 

Commercial  Make 

MIPS  R3010 

Weitek  3364 

TI  8847 

KVKWiiiL^VTail 

Cycle  Time(ns) 

40 

50 

30 

20 

Cycles/add 

2 

2 

2 

2 

Cycles/mult 

5  : 

2 

3 

2 

Cycles/divide 

19 

17 

11 

3 

Cycles/sq  root 

- 

30 

14 

? 

Data  on  MIPS  R3010,  Weitek  3364,  and  TI  8847  was  obtained  from  the  computer 
architecture  book  by  Hennesey  and  Patetrson  [Ref.  30]  while  the  data  for  AT&T 
DSP32c  is  obtained  from  DSP32c  data  manual  [Ref.  31].  Fig  5.1.  gives  a  schematic 
of  the  2nd  order  HR  filter  hardware  scheme.  It  has  a  coefficient  memory  which  can 
be  set  by  an  external  device.  Close  examination  reveals  that  Multiplications  A  can 
be  performed  simultaneously  while  Multiplication  B  and  other  additions  are  to  be  in 
sequence.  Fig  5.2.  is  the  basic  building  block  for  the  adaptive  filter  and  is  identical  for 
all  the  sections.  This  is  offers  an  advantage  in  hardware  complexity  over  the  Kwan  & 

Martin  approach.  Filters  1  and  2  have  transfer  functions  of  Hbp{z)  while  filter  3  has  a 
Cx  do*CL  • 
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transfer  function  Hgt{z).  Throughput  rate  is  primarily  determined  by  the  computing 
time  of  the  2nd  order  HR  filter.  The  overall  block  diagram  of  the  adaptive  filter  for  a 
3-Notch  cancellation  was  given  in  Fig  5.3.  This  architecture  remains  same  either  for 
a  Floating  point  or  Fixed  Point  arithmetic. 

1.  Time  Budget 

In  the  architecture  of  Fig  5.1,  5.2,  and  5.3  the  most  vital  elements  are  the 
UR  filter,  Controller  and  the  Forgetting  filter.  The  HR  filter  corresponds  to  Hbp{z). 

#■ 

The  controller  has  to  update  the  par2uneter  via  equation  (3.17)  which  calls  for  2 
*  multiplications,  1  addition,  and  1  division.  Sirnilairly  the  Forgetting  Filter  has  to 

implement  equation  (3.19)  calling  for  2  multiplications  and  1  addition.  Results  are 
tabulated  as  given  below 


TABLE  5.2 


Element 

of  Effective 
Multiplications 

#  of  Effective 
Additions 

#of 

Divisions 

HR  filter 

2 

3 

nil 

Controller 

2 

1 

1 

Fogetting 

Filter 

2 

1 

nil 

The  maximum  computing  time  is  at  the  HR  filter.  Using  the  AT&:T  DSP32c  proces¬ 
sor,  it  takes  5x2  «=  10  cylces  (ie  200ns)  for  implementation  the  desired  HR  filter.  This 
corresponds  to  a  Throughput  rate  of  bMhz.  This  is  the  best  that  could  be  achieved 
with  the  existing  commercial  floating  point  processors.  However  processors  like  the 
HDSP661 10  of  Honey  well  make  claim  an  ALU  speed  of  10ns  with  block  floating  point 
operations  can  conveniently  implement  at  \0Mhz  throughput  rate,  without  pipelin¬ 
ing  the  HR  filter.  With  pipelining  the  throughput  rate  would  increase  dramatically, 
but  will  never  exceed  the  maximum  throughput  rate  of  lOOAffi;. 
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2.  VLSI  Approach 

Taking  advantage  of  the  fact  that  the  entire  hardware  can  be  generated  by 
a  simple  repetitive  use  of  the  building  block,  strong  consideration  should  be  given 
to  designing  a  VLSI  chip  for  Fig  5.2  rather  them  implementing  the  hardware  via  the 
existing  commercial  chips.  The  main  rationale  behind  this  idea  is  that  reliability 
of  the  system  is  very  important  in  EW  equipment.  Also  VLSI  has  the  potential  of 
ultimate-low  cost. 
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VI.  CONCLUSIONS 


In  this  thesis  the  problem  of  suppressing  narrow-band  interference  was  ad¬ 
dressed.  The  problem  specific  to  the  Electronic  warfare  scenario  was  kept  in  mind 
while  solving  the  problem.  The  new  algorithm  [Ref.  26]  was  derived  and  simulated. 
This  new  algorithm  was  tested  against  various  signals  amd  signal  conditions.  The 
results  are  highly  encouraging. 

Some  aspects  of  pipelining  were  discussed.  Reduced  hardware  complexity  of  the 
New  Algorithm  is  an  important  advantage  over  the  ewlier  algorithm  by  Kwan  and 
Martin  [Ref.  8].  A  possible  hardware  scheme  for  implementing  the  new  algorithm 
was  discussed.  It  was  observed  that  with  the  existing  commercially  available  floating 
point  processors  a  throughput  rate  of  hMhz  is  achievable  while  using  processors  like 
HDSP66110  with  an  ALU  speed  of  10ns  it  is  possible  to  achieve  a  throughput  rate 
of  IQMhz. 

A.  FUTURE  WORK 

Future  directions  of  work  are 

i)  VLSI  design  of  the  system  with  an  architecture 
as  indicated  earlier  or  a  similar  architecture. 

ii)  Pipelining  of  2nd  Order  HR  filter  suited  for 
this  application 

iii)  Design  using  recursive  DFT 

The  above  areas  of  work  are  the  logical  extensions  of  this  work.  However  a  radically 
different  approach  for  the  same  problem  using  adaptive  filtering  in  frequency  domain 
[Ref.  10]  should  also  be  considered. 
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APPENDIX  A 


dimension  faray(5000) 
dimension  uaray(5000) ,yaray(5000) 
dimension  fhat (5000) ,ph( 5000) 
dimension  a(90),b(90) 

open  (unit=9 , f ile= ' spk. dat ' , status® • new ' ) 
c 

c  This  programme  generates  input  output 
c  sequence  by  exciting  a  linear  system 
c  defined  by  the  numerator  polynomial 
c  B(z)  and  denominator  polynomial  A(z) 
c  this  data  is  stored  in  uarray  and  yarray 
c  Subsequently  an  ARMA  (pole, zero)  is  used 
c  to  fit  this  data, 
c 

n=1024 

ix®l 

yk=rand(ix) 

ix=0 

kk=8 

c 

c  generate  sequence  uk  and  yk 
c 

ip»3 

iz«2 

a(l)®1.7 
a(2)«-1.53 
a(3)»0.648 
b(l)*l,0 
b(2)-0,6 
do  10  k>l,n 
call  rbn(uk,ix,k) 
call  sy8tem(a,b,uk,yk,ip,iz,k) 
uaray(k) ^uk 
yaray(k)«yk 
10  continue 

save  the  sequence  in  the  arrays  uaray  and  yaray 

call  8pktrm(yaray,faray) 
c 

c  spectrum  of  the  output  sequence  yaray  is  computed 
c  using  FFT  programme  and  stored  in  faray  for 
c  comparision 
c 

ip«4 

iz«l 

c 

c  fit  an  AHHA(ip, iz)  model  for  the  given  data 
c  this  is  an  ordinary  least  squares  algorithm 
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c 

print  * ,  ’  inpur-ip-iz“>  • 
read  *,ip,iz 
do  20  k-l,n 
uk»«aray (k) 
yk=yaray(k) 

call  araa(a,b,uk,yk,lp,iz,kk,k) 
c  print  'nuin->' ,  (b(i)  ,i»l,iz) 

c  print  *den->  * ,  (v^(i)  ,i«l,ip) 

20  continue 

print  *, 'nuin->’ ,  (b(i) ,  i«l,iz) 
print  •den->' ,  (aH)  f  i*l»ip) 
c 

c  Sampling  of  the  Z-transform.  given 

c  coefficients  a  . ^ .  and  b  . 

c  corresponding  H(z)  *>  B(z)/A(z) 
c  is  evaluated  on  the  unit  circle 
c  for  obtaing  U>e  parametric  specxrum 
c 

call  zsrople(a*b,fhat,ph,ip«iz) 
c 
c 

do  30  i»l,n/2 

c  print  *,i,  >',fhat(i)  ,faray(i) 
write  (9,*)  i,fhat(i) ,fardy(i) 

30  continue 
stop 

c 

c 

subroutine  zsmple(aib«r,ph, ip,iz) 

dimension  a(90),b(90) 

dimension  r(2046)fPh(204d) 

complex  z » ui , f n « f d , omega , delw , spec 

pi**atan(1.0) 

pi«4.0*pi 

ui«{0. 0,1.0) 

delw* (pi/1024 . 0) *ui*2 . 0 

omega*(0. 0,0.0} 

do  100  k«l,1024 

z«exp( -omega) 

fd-‘{0.0,0.0) 

do  10  i«l,ip 

10  fd*fd+a(ii*(za*(-i)) 
fd-l.O-fd 
fn>»(0. 0,0.0) 
do  20  i*l,iz 

20  £n*fn^b(i)a(z«*(-i)) 

8pec«>fn/fd 

r(k)«ab5(spec) 

S|^cr*real  (  spec) 
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speci^aimag ( spec) 
ratio»speci/specr 
ph (k) =atand (ratio) 
omega«omega+delw 
c  print  * , '  — z-oin-> '  ,z, omega 
100  continue 
return 
end 
c 
c 

function  atand(x) 
piBatan(l.O) 
pi=4 . 0*pi 

c  xrad«atan(x) 

atandsxrad* ( 180 . o/pi ) 
return 

*  end 

c 

c 

subroutine  rbn(b«ix,k) 

z«rand(ix)-0.5 

if(z.ge.O.O)  b=1.0 

if(z.it.O.O)  b»-1.0 

return 

end 

c 

o 

subroutine  6ystem(a«b,uk«yk,lp,i2«k) 
dimension 

«  zkar(90) ,zkma(90} , 

4  a(90),b(90) 

if  (k.ne.l)  go  to  250 
do  230  i**l,ip 
230  zkar(i)"0«0 

do  240  i>l»iz 
240  zktta(l)«‘0.0 

250  continue 

do  30  i*ip,2,-l 

*  30  zkar(i)«zkar(i-l) 

zkarU)*yk 

,  do  31  i»iz,2,-l 

31  zksia(i)>«zkina(i<^l) 
zkDa(l)^uk 
c 
c 

yk**0 .0 
do  200 

200  yk«yk4a(i)*2kar(i) 
do  210  {•‘Itiz 
210  yk“yk+zk®a(i)*b(i) 
return 
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end 

c 

c 

subroutine  delay {uX,yk,idelay,3c) 
dimension  d(500} 
if  (k.ne.l)  goto  10 
do  20  i=l,500 
20  d(i)»0.0 

10  continue 

do  30  i»500,2»-l 
30  d(i)»=d(i-'l) 

d(l)»uk 
yk®d(idelay+l) 
return 
end 
c 
c 

subroutine  spktrm{taray,faray) 
complex  u(SOOO) 

dimension  taray(5000) ,faray(5000) 

m*10 

n*1024 

rn«sn 

do  20  i«*l,n 
u(i)*»taray(i) 

20  continue 

call  pf{t(u,B,n) 
do  30  i“l,n/2 
f aray ( 1 ) «real (u ( i) ) 

30  continue 
return 
end 
c 
c 

o  *«***************«I*  *<?#•**<*«******* 

c  * 

c  Subroutine  for  computing  DFt  of  * 

c  an  array  *a*  is  complex  and  a  pair 
c  cf  numbers  are  to  be  specified  * 

c  for  each  point  * 

c  M  is  the  2  power  index  * 

c  say  ID»10  then  number  of  points  * 

c  are  1024  * 

c  after  computation  fft  is  kept  * 

c  in  the  same  complex  array  *a*  • 

c  e 

C  ********************************** 

c 

subroutine  pfftCa^B^n) 
complex  a(5000) «UtV,t 
complex  uifUr 
pi«3 . 141592653509793 
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n=2**iii 

nv2=n/2 

mnl=n-l 

j=l 

do  10  i=l,mnl 
if(i.ge.j)  goto  20 
t»a(j) 
a(j)=a(i) 
a(i)«t 
20  k-nv2 

30  if(k.ge.j)  goto  10 
j=j-k 
k*k/2 
goto  30 
10  j*j+k 

do  40  1=1, m 
le=2**l 
lel=le/2 
u=(l. 0,0.0) 

w=cinplx(cos(pi/lel)  ,-sin(pi/lel) ) 
do  40  j=l,lel 
do  50  i=j,n,le 
ip=i+lel 
t=a(ip) *u 
a(ip)=a(i)-t 
50  a(i)=a(i)+t 

40  u=u*w 
rn=n/2 

ui«(0. 0,1.0) 
ur=(l.0,0.0) 
do  60  i=l,n/2 
a(i)=a(i)/rn 
a(i+512)=a(i+512)/rn 
areal=real(a(i) ) 
aiingn=aiiaag(a(i) ) 
anagn=abs(a(i) ) 
aphase=atand (alngn/areal) 
a  ( i)  =ur*ainagn+ui*aphase 
areal=real (a ( i+512 ) ) 
aimgn=aiiDag(a(i+512) ) 
ainagn=abs  ( a  ( 1+512 ) ) 
apha  se=a tand ( a imgn/ area 1 ) 
a ( i+512 ) =ur*amagn+ui*aphase 
60  continue 
return 
end 
c 

c  $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

c 

c 

subroutine  anaa(a,b,uk,yk, lp,iz,kk,k) 
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dimension 

&  zkar(90)  ,z)aaa(90) , 

&  zk(90) ,a(90) ,b(90) ,c(90) 

if  (k.ne.l)  go  to  250 
n=ip+iz 
ipp=ip+l 
do  230  i=l,ip 
230  a(i)=0.1 

do  240  i«l,iz 
240  b(i)=0.1 

250  continue 

do  30  i«ip,2,-l 

30  zkar(i)=zkar(i~l) 
zkar(l)=ykold 
ykold*yk 

do  31  i=iz,2,-l 

31  zkma(i)=zkma(i-l) 
zkma(l)=uk 

do  32  i=l,ip 
c(i)«a(i) 

32  zk(i)»zkar(i) 
do  33  i=l,iz 
c(i+ip)«b(i) 

33  zk(i+ip)«zkma(i) 

call  syseq(c,yk, zk,n,kk,k) 
do  34  i=l,ip 

34  a(i)=c(i) 

do  35  i«l,iz 

35  b(i)=c(i+ip) 

return 

end 

c 

c 

c 

subroutine  syseg(a,yk,zk,n,kk,k) 
c 

c  solves  system  of  equations  with 

c  I  of  equations  more  than  unknowns 

c  using  linear  reggression. 

c  t 

c  equation  is  yk=zk  *  a 

c  where  'a'  is  n  vector  to  be  estimated, 

c 

C  *4t4t***A*4rAj(Aik4i****Ait4t4f******AA********** 

C 

dimension 

6  del(90) ,phat(90,90) ,zk(90) ,y(90) ,a(90) 

if  (k.ne.l)  go  to  250 
alpha=1.0 
250  continue 

call  lnv(phat,zk, alpha, n,kk,k} 
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* 


t 


ykhat=0 . 0 
do  40  i=l,n 

40  ykhat=ykhat+zk(i)*a(i) 
er=yk-ykhat 
do  60  i=l,n 
60  y(i)=er*zk(i) 
do  80  i=l,n 
del (i) =0.0 
do  80  j=l,n 

80  del(i)=phat(i, j)*y(j)+del(i) 
if  (inod(k,kk)  .ne.O)  goto  70 
c  print  >',k 

c  print  100,  ( (phat(i, j) ,i=l,n) , j=l,n) 

c  print  130,  (del(i) , i=l,n) 

c  print  160, er 

c  print  170,  (zk( j ) , j=l,n) 

70  continue 

do  50  i=l,n 
50  a(i)=a(i)+del(i) 

100  format(2x, 'phat' ,2x,/,4el6.6) 

120  format(2x, 'zk' ,2x,/,4el6.6) 

130  forniat(2x, 'del*  ,2x,/,4el6.6) 

160  format (2x, 'ekl ', 2x,el6. 6,/) 

170  format(2x, 'zk' ,2x,/,4el6.6) 
return 
end 
c 
c 

subroutine  inv ( phat , zk , alpha , n , kk , k) 

dimension  phat(90,90) ,zk(90) ,q(90) ,r(90) ,delp(90,90) 

if(k.ne.l)  goto  5 

do  6  i=l,n 

do  6  j=l,n 

if(i.eq.j)  phat(i, j)=1.0 
if(i.ne.j)  phat (i,j) =0.0 
6  continue 

5  continue 

do  90  i=l,n 

do  90  j=l,n 

phat ( i , j ) =phat ( i , j ) /alpha 
90  continue 

do  10  i=l,n 

q(i)«0.0 
do  10  j=l,n 

10  q(i)=phat(i, j)*zk(j)+q(i) 

8uin=1.0 

do  20  i=l,n 

20  suin=sum+zk(i)  *q(i) 
do  30  j=l,n 
r(j)=0.0 
do  30  i=l,n 

30  r(j)=phat(i, j)*zk(i)+r(j) 
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do  40  j=l,R 
do  40  i=l,n 

40  delp(i, j)=q(i)*r(j) 
do  50  j=l,n 
do  50  i=l,n 

50  phat(i,j)»phat(io)-(delp(i,j)/sum) 
return 
end 
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c 

c 

c 


APPENDIX  B 


simulation  of  michal  paper  IEEE 

dimension  alcl(lO)  ,ak2(10)  ,s(10)  ,xout(10)  »theta(10) 
dimension  a)cs(lO)  ,wp(10)  ,ak3  (10) 
dimension  array (4, 4000) 

open (  unit=9 , file= ’martin. dat ' , status= ' new ' ) 
open(  unit=8,file=' mart. dat' , status®' new ' ) 

print  *, '-input-SNR-in-dB - >' 

read  *,snr 
c  if  (snr.ne.O)  goto  200 
c  var=o . 5 
c  goto  210 

var=1.0/(2.0*(10.0**(snr/10.0) ) ) 

210  continue 

print  * , '  — variance— > ' ,  var 
n=l 

print  * , ' -input-nh — > ' 
read  *,n 

print  *, '-input-# -of -waves — >' 
read  *,nd 

print  * , '  -input-angles— > ' 
read  (theta(i) , i=l,n) 
r=0.9 

pi=4.0*atan(1.0) 
rad»pi/180.0 
amu»0 . 0 

c  do  50  kl"l,50 
amu»0.015 

print  * , '  -input-step-size—  0 . 015— > ' ,  amu 
read  *,amu 

print  *, ' — angles — >' , (theta (i) ,i»l,n) 

print  *, ' — input — l-for-step-chango— else-O — >' 

read  *,ic 

if  (ic.eq.l)  print  input-step-in-degrees— >' 

if  (ic.eq.l)  read  *,step 

do  40  i®l,n 

ak2(i)-(l-r*r) 

ux>theta (i) *rad 

value«cos (ux) 

c  akl  ( i )  -Kll  *8qrt  ( 1 . 0-f  r*r-2  *  0*r*value) 

akl ( i ) -sqrt ( 1 . O+r *r ) 
aks  ( i)  >sqrt  ( 1 . 0’fr^r-a .  o*r*value) 

40  continue 
o  amu-anu-f 0 . 0001 

do  10  k«l,6500 

call  system (yk, wk, theta , step, var ^ ic , nd , k) 
call  f liters (akl, ak2, aks, 8, en,xout,yk,n,k) 
call  increment (akl ,ak2 , aks, s, en, amu, avg, n, k) 
do  30  i*l,n 
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xx=2.0*sqrt(l-0.5*ak2 (i) ) 
wp(i)=2.0*asin(akl(i)/xx) 
c  wp(i)*wp(i)/rad 
wp(i)«cos(wp(i) ) 

30  continue 

do  60  i=l,n 

ak3 ( i) =abs (akl ( i) ) 

60  continue 

kk=mod(k,100) 

ll=mod(k,10) 

if(kk.eg.O)  print  ' — akl — >' , (ak3 ( j) , j-l,n) 

if(ll.eq.O)  write(9,*)  k, (ak3 (j) , j=l/n) 
if(kk.eq.O)  print  *,• — aks — >' , (aks(j) , j=l,n) 
c  if(kk.eq.O)  print  *,k, '-avg->' ,avg 

c  print  (wp(i) ,i=l,n) 

c  write(9,*)  k, (wp(i) ,i=l,n) 

if (k.gt. 1500. and. k. It. 2500)  write(8,*}  k,yk,en 
if (k.gt. 5000. and. k. It. 6000)  write(8,*)  k,yk,en 
10  continue 

c  print  *,  kl,amu,avg 

c  write (9,*)  kl,amu,avg 

50  continue 

close (9) 
stop 
end 
c 
c 

subroutine  increment ( akl , ak2 , aks , s , en , amu , avg , n , k) 
dimension  akl (10) ,ak2(10) , aks (10) ,s(10) ,ss(10) 
dimension  dec (10) , flag (10) , fading (10) 
if  (k.ne.l)  goto  30 
c  amu>»0.001 

avg=0.0 
fading(l)=0.9 
fading (2) =0.9 
fading (3) =0.9 
fading (4) =0.9 
30  continue 

do  10  i=l,n 

ss (i) “fading (i) *ss(i)+s(i)*s(i) *(1.0-fading(i) ) 

10  continue 

c  do  31  i=l,n 

c31  print  *, '-en-s-ss->' ,en,s(i) ,ss(i) 

c  sum=0.o 

do  20  i=l,n 

dec  (i)=amu*en<'s(i)/(ss(i) +0.0001) 
akl(i)=akl(i)-dec(i) 

c  sum=sum+(abs(aks(i)-akl(i) ) )/(aks(i) ) 
call  stability (akl, ak2,f lag, n) 
akl ( i ) «akl ( i ) +f lag ( i ) *dec ( i ) 

20  continue 

realk=k 
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r 


» 


# 


k 


avg=avg*  ( ( realk-1 . 0 )  /realk)  +suin/realk 

return 

end 


svibrout Ine  filters  ( akl ,  ak2 ,  aks ,  s ,  en ,  xout ,  xin ,  n ,  kl ) 
distension  w(10|3)  ,a(10,2)  ,b(10,2)  «e(ll)  ,xbp(10} 
distension  g(10,3)  ,gbp(10)  ,u(10) 
distension  q(10,3)  ,s(10) , xout (10)  »c(10) 
distension  akl(lO)  ,ak2(10)  ,aks(10)  ,as(10) 
c 

if  (ki.ne.l)  goto  200 

c  open (  unit=9 ,  f  ile= ' startin. dat ' , status® ' new ' ) 

200  continue 
i«n 

do  40  i=l,n 

a(i,l)=-(akl(i)*akl(i)+ak2(i)-2.0) 

as(i)B-(aks(i)*aks(i)-t-ak2(i)*‘2.0) 

a(i,2)=-(1.0~ak2(i)) 

b(i,l)=-0.5*ak2(i) 

b(i,2)=0.5*ak2(i) 

c(i)=-akl(i)*ak2 (i) 

40  continue 
c  do  41  i®l,n 

c  print  *, '-kl-k2->* ,akl(i) ,ak2(i) 

c4l  print  *, ’-a-b-c->' ,a(i,l) ,a(i#2) ,b(i,l) ,b(i,2) ,c(i) 


e(l)*xin 

i«n 

do  10  j=l,n 

w(i,3)*w(i,2) 

w(i;2)'=w(i,l) 

w(i,l)*a(i,l)*w(i,2)+a(i,2)*w(i,3)+e(j) 

xbp(j)-b(i,l)*w(i,l)+b(i,2)*w(i,3) 

e(j+l)«xbp(j)+e(j) 

i»i-l 

10  continue 
en®e{n+l) 


i®n 

do  20  k®l,n 
u(k)»en-xbp(k) 
g(i,3)-g(i,2) 
g(i,2)-g(i,l) 

g(i,l)-a(i,l)*g(i,2)+a(i,2)*g(i,3)-fu(k) 

gbp(k)-b(i,l)*g(i,l)+b{i,2)*g(i,3) 

q(k,3)«q(k,2) 

q(k,2)®q(k,l) 

q(k,l)«a(i,l)*q(k,2)+a(i,2)*q(k,3)+gbp(k) 
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20 

c 

c 

c 

c 

c 


c 

c 

c 


100 

c 

c 


10 

o 

c 


c 

c 


20 

10 


30 


s(i)=c(i)*q(k,2) 

i=i-l 


continue 

kk=inod(ki,10) 

if(3ck.eq.0)  write(9,*)  ki,en,a(l,l) /as(l)  ,a(2,l) 
if(kk.eq.O)  print  *,  ki, (a(i, 1} , i=l»n) 

if(kk.eq.O)  print*,  ki, (as(i) ,i=l,n) 

if (ki.ge.1800)  write(9,*)  ki,e(l) ,en,xbp(l) ,a(l, 
if (ki.ge.1800)  print  *,  ki,e(l) ,en,xbp(l) ,a(l, 
return 


end 


subroutine  sy steia  (yk ,  wk ,  theta ,  step ,  var ,  ic ,  n ,  k) 
dimension  theta (10) ,yout(l0) ,yz(10) 
if  (k.ne.l)  goto  100 
continue 

call  bpf (theta, yz.n,k) 

call  waves ( theta, yout,n,k) 

call  file(yk,theta,step,ic,n,k) 

sum^O.O 

wk*0 . 0 

do  10  i«l,n 

6uoi[*6um*yout  ( i ) 

wk'-wk+yz  ( i) 

continue 

call  gnoise(g,var,k) 

yk»8un*g 

yk»wk+g 

ykeyk+g 

return 

end 


subroutine  waves (theta, yout,n,k) 

dimension  zkar(10,3) ,theta(10) ,yout(10) ,al(10) 

if (k.ne.l)  goto  10 

pi«(atan(1.0) )*4.0 

rad“pi/180.0 

do  20  i««l,n 

2kar(i,l)*0.0 

2kar (i , 2) *“8in (theta ( i) *rad) 

al ( i ) -2 . 0*cos (theta ( i ) *rad) 

continue 

continue 

do  30  i«l,n 


2kar(i,3)«2kar(i,2) 
2kar(i,2)>zkar(i,  1} 
*kar(i,l)-al(i>*zkar(i,2)-zkar(i,3) 
yout (i)«2kar(i,l) 
continue 


,as(2) 

1) »as(l) 
1) ,as(l) 
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» 


# 


i 


return 

end 

c 

c 

subroutine  gnoise(g,s,lc) 
if()c.ne.l)  goto  10 
ix=l 

yfl’=rand(ix) 

ix*0 

10  continue 
6um=0 . 0 
do  20  i»l,12 
yfl*rand(ix) 

20  suffiasurn't^yfl 

gs(SU]n->6.0)*6 

return 

end 

c 

c 

subroutine  stability (akl,ak2, flag, n) 
dimension  akl(lO) ,ak2(10) , flag (10) ,a(10,2} 
do  30  i>l,n 
flag(i)-0.0 

a(i,l)  — (akl(i)*akl(i)+ak2(i)-2.0) 
a(i,2)»-(l,o-ak2(i)) 
des»a(i,l)M(i,l)+4.0*a(i,2) 
if  (des.lt. 0.0)  goto  10 
desoO . 5*6grt (des) 
root l«0 . 5*a ( 1 , 1 ) +des 
rl«ab6(rootl) 
root2«0 . 5*a ( i , 1)  -^es 
r2'«ab6(root2) 
if  (rl.ge.l.O)  goto  20 
if  (r2.ge.l.O)  goto  20 
c  print  * ,  •  *“real -roots— >  * ,  rootl ,  root2 
goto  30 

20  flag(i)»>1.0 

goto  30 
10  deso-des 

y«‘0.5*sgrt(des) 

x-0.5*a(i,l) 

xx-'ab8(x) 

if  (xx.le.l.Oe-6)  x-0. 0000001 
angle«(atan(y/x) )*57.3 
radiusesgrt ( x*x4y •y ) 

c  print  • ,  ’  —complex— roots— >  * ,  radius ,  angle 
if (radius. ge. 1.0)  goto  20 
30  continue 
return 
end 


c 

c 
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c 

c 

c 


$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 


subroutine  bpf( theta, yout,n,k) 

dimension  zkar(10,3) ,a(10,3} ,zkma(10,4) ,b(10,3) 
dimension  theta(10} ,g(10) ,yin(10) ,yout(10) ,uk(10} 
if  (k.ne.l)  go  to  10 
pi«4.0*?itan(1.0) 
rad*pi/  180.0 
do  15  i=l,n 
zkar(i,l)a0.0 
15  zkar(i,2)«0.0 

r«0.985 

c  print  * , ' — input — r — > • , r 
c  read  * , r 

do  25  i-l,n 

a(i,3)-r*r 

g(i)»0.5*(i-r*r) 

b(i,l)«g(i) 

b(i,3)«-g(i) 

a ( i , 2 ) «2*r*cos (theta ( i) *rad) 

25  b(i,2)-0.0 

10  continue 
8ave«var 
do  35  i>l,n 
var«2 . 0 

call  9noise(gk,var,k) 
yin(i)«gk 

zkar ( 1,3) •zkar (1,2) 
zkar(l,2)«zkarU,l} 
zkna ( 1 f  3 ) ozkma (1,2) 
zkBa(l,2)-zkma(i,l) 
zkffla(l,l)«yln(l) 

uk(l)«z)cna(i,  1)  *b(i,  l)~zkaa(l,2)*b(l,2)’fzkffla(l,3)  *b(l,3) 
zkar(l,l)«zkar(l,2)*a(l,2)‘‘zkar(l,3)*a(l,3)4>uk(l} 
c  print  *, '“-f  11— >',af  11, theta, r 
c  print  *, '—f  11— >' ,zkar 
35  yout(l)ozkar(l,l) 
var-save 
return 
end 

c  $$$$S5$$S$$$$$$S$$$$$$$$$$$S$$$$$$$$$$$ 
o 

subroutine  file (yk, theta, step, tc,n,k) 
dimension  output ( 8000) , theta (10) 

If  (k.ne.l)  goto  10 

call  gendata ( output , theta , step , Ic , n , k) 

10  continue 

yk>‘OUtput(k) 

return 

end 
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c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


f 


9 


i 


c 

subroutine  gendata (output , theta , step , ic , nn , k) 

DIMENSION  YSAMP(65536) ,data(4,8000) , output (8000) 
dimension  theta (10) 

c  OPEN  (UNIT=9 , FILE*  * bpsk.dat • , STATUS* • NEW • ) 
n»7000 
amag*l.4l4 
tcarr*4 . 8 
tchip»7 . 0 
ichip*0 
tdata*n 
tdelay*o.o 
1*9999 
tchip*200.0 
do  10  i*l,nn 
tcarr*3  60 . o/theta ( i ) 

call  bp8k(n,amag»tcarr,tchip,ichip,tdata,tdelay,l,ysamp) 

do  20  j«l,n 

if  (io.eq.O)  goto  80 

if  U.eq.2)  goto  SO 

data(i,  j)-y8ainp(j) 

goto  60 

50  if  (j.gt.3000)  y6amp(j)*0.0 

80  dataU#  j)'‘ysaap(j) 

60  continue 
20  continue 

print  *,i, ,theta(i) ,tcarr 
10  continue 

tcarr«360. 0/90.0 
tchip«8.0 

call  bpsk (n , amag , tcarr , tchip, ichip, tdata , tdelay , 1 , ysamp) 

do  30  i*l,n 

x*0.0 

do  40  j»l«nn 
40  it“)t+data(j,i) 

output  ( 1 )  *x-t-y  samp  ( i ) 
c  vrite(9,*)  X 
30  continue 

if  (ic.eg.O)  return 
nl*n-3000 
tchip*200.0 
thetal*theta  (2 )  -t^step 
tcarr*360 . O/thetal 

cal 1  bpsk ( nl , amag , tcarr , tchip, ichip , tdata , tdelay , 1 « ysamp) 
j“l 

do  70  i*3000,n 
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output ( i ) "Output ( i ) +y samp ( j ) 

j«j+l 

continue 

close (9) 

return 

end 


this  programme  was  modified  to  suit  the  martin 

programme 

subroutine  bpsk 

&  (n,  asiag ,  tcarr ,  tchip«  ichip ,  tdata ,  tdelay ,  11 ,  ysamp) 


c 

C  THIS  PROGRAM  GENERATES  SAMPXiES  OF  A  DIRECT  SEQUENCE  &I**PHASE 
SHIFT 

C  KEYED  SPREAD  SPECTRUM  SIGNAL.  THE  '’INFORMATION"  BITS  AND  THE 
C  SPREADING  SEQUENCE  USED  ARE  RANDOMLY  GENERATED.  PARAMETERS 
REQUIRED 

C  FOR  OPERATION  ARB: 

C  N  "  NUMBER  OF  SAMPLES  GENERATED 

C  FOR  CONSISTENCY  WITH  USE  BY  AN  FFT 

C  ALGORITHM,  N  SHOULD  BE  AN  INTEGER 

C  POWER  OP  2  —  TYPICALLY  1024 

C  NOTE:  IF  N>1024  DIM  OF  YSAMP  MUST  CHANGE 

C  MAG  •  MAGNITUDE  OF  CARRIER  WAVEFORM 

C  TSDELAY  «  NUMBER  OF  SAMPLE  TIMES  DELAYED  FROM 

C  t  «  0  BEFORE  BEGINNING  SAMPLES 

C  AN  ARBITRARY  VALUE  THAT  ALLOWS  SOME 

C  FLEXIBILITY  IN  SAMPLING.  SHOUU)  BE 

C  ZERO  FOR  SAMPLING  AT  t-0. 

C  TDATA  "  NUMBER  OP  SAMPLE  TIMES  IK  ONE  DATA  BIT 

C  TCUXP  "  HUMBER  OF  SAMPLE  TINES  IN  ONE  BIT  OF 

C  SPREADING  CODE 

C  TCARR  -  NUMBER  OF  SAMPLE  TIKES  IN  ONE  CYCLE  OF 

C  CARRIER  FREQUENCY 

C  T5AHP  "  DURATION  OF  A  SAMPLE  TIME.  IN  GENERAL 

C  TSAMP  WILL  ALWAYS  BE  "  1.0,  SINCE 

C  TIKE  SAMPLED  VALUES  BECOME  STATIC 

C  WHEN  STORED  AND  CAN  BE  SCALED  LATER. 

C 

C  THE  ALGORITHM  USED  TO  GENERATE  THE  DS-BPSK  SIGNAL  RELIES  UPON 


TDATA 

TCHIP 

TCARR 

TSAMP 


BUILTIN  FORTRAN  RANDOM  NUMBER  GENERATOR  RAND(L)  TO  PRODUCE  THE 


C  INFORMATION  BIT  STREAM  AND  THE  SPREADING  SEQUENCE  CODE.  AT 
THE  TIME 

C  OF  EACH  SAMPLE,  THESE  TWO  BITS  ARE  MODULO  2  ADDED  TO  PRODUCE 
A  CHIP 

C  BIT.  THIS  BIT  IS  THEN  USED  TO  SHIFT  THE  PHASE  OF  THE  CARRIER 
BY 
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C  180  DEGREES  IF  1  OR  BY  0  DEGREES  IF  0.  THE  RESULTANT  VALUE  OF 
THE 

C  COSINE  FUNCTION  IS  MULTIPLIED  BY  THE  AMPLITUDE  OF  THE 
WAVEFORM,  THEN 

C  STORED  IN  FILE  BT.DAT  IN  THIS  FILE  DIRECTORY.  THE  PROCESS  IS 
CONTIN- 

C  UED  UNTIL  N  SAMPLES  ARE  GENERATED.  NOTE  THAT  FIRST  LINE  OF 
BT.DAT 

C  CONTAINS  N,  DSBPSK,  DATE  AND  TIME  AS  ID  FIELD  FOR  FILE. 

C 

C  W.R. TUCKER  9-MAY-83 

C  CONVERTED  PROGRAM  TO  ALLOW  FOR  M-SEQUENCE  CHIP 

C  M-SEQUENCE  GENERATOR  REQUIRES  DATA  FILE  FSR.DAT  WITH 

(  C  PARAMETERS  FOR  FEEDBACK  SHIFT  REGISTER. 

C  W.R. TUCKER  4-AUG-83 

C 

»  C  INITIALIZE  —  I  -  INFO  BIT  NUMBER,  J  -  SP  SEQ  CODE  BIT  NUMBER 

C  L  IS  ARGUMENT  FOR  RAND  (L)— A  DIFFERENT  SEQUENCE  MAY  BE 
C  GENERATED  BY  INITIALIZING  L  WITH  DIFFERENT  VALUES. 

C 

C  Modified  by  HHL  for  installation  on  ULTRIX  ECE  VAX,  12/90 
c  Writes  output  in  file  bpsk.dat 
c 

DIMENSION  YSAMP(65B36) 

INTEGER  I,J,L,INFO,CHIP,lFSR,IFB0(100) ,IVAL(100) 

REAL  TDATA,TCHIP,TCARR,TDELAY,TSAMP,KAG,yARG,YSAMP,PI 
PI  »  3.14159265 
TSAMP  *  1.0 
I  «  0 
J  -  0 
L  «•  0 

WRITE (6, 900) 

900  FORMAT (*  GENERATION  OF  DIRECT  SEQUENCE  BPSX  SIGNAL*, 

A  *  IN  FILE  bpsk.dat*} 

WRITE (6, 1000) 

1000  FORMAT (*  I  SAMPLES  TO  GENERATE  •  *,$) 

C  R£AD(S,  *)N 
print  ♦,n 
WRITE (6, 1001) 

^  1001  FORMAT (♦  MAX  AMPLITUDE  OF  SAMPIil  VALUES  (R)=  *,S) 

C  READ(S,  *}  MAG 

^  ttag«aitiag 

print  *,Ba9 
WRITE (6, 1002} 

1002  FORMAT (*  I  SAMPLES  PER  CARRIER  CYCLE  (R)«  *,$} 

C  R£AD(5,  *}  TCARR 

print  *,tcarr 
WRITE (6, 1003) 

1003  FORMAT(*  I  SAMPLES  PER  CHIP  BIT  (R)  *>  *,$) 

C  R£AD(S,  *)  TCHIP 

print  *,tchip 
WRITE (6, 1020) 
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1020  FORMAT(’  ENTER  (0) : RANDOM  CHIP,  OR  (1) ; REPEAT  M-SEQ  CHIP 

C  '  READ(5,  *)  ICHIP 
print  *,ichip 
WRITE (6, 1004) 

1004  FORMAT ('  #  SAMPLES  PER  INFO  BIT  (R)  =  ',$) 

C  READ(5,  *)  TDATA 

print  *,tdata 
WRITE(6,1005) 

1005  FORMAT ('  «  DELAYS  BEFORE  SAMPLING  (R)-  ',$) 

C  READ(5,  *)TDELAY 

print  *,tdelay 
WRITE (6, 1006) 

1006  FORMAT(‘  RANDOM  NUMBER  SEED  (14) - ■>  ’,$) 

1»11 

C  R£AD(5,  *)L 
print  *,1 

c  Initialize  RAND  by  emailing  with  input  non'*zero  L. 
o  Subsequent  calls  will  be  with  L  »  0. 

RANDOM»RANi:^(L) 
c  debug 

C  WRITE{6,*)  *  SEED  A)m  RANDOM  NUMBER  RETURNED: ',L, RANDOM 
1^0 

WRITE<6,i010) 

1010  FORMAT (/IX.  ’ - WORKING-— - ^-« ) 

IF  (ICHIP  .EQ.  0)  GO  TO  S 

OPEN  (UNIT  «  1,  FILE  -  ♦FSR.DAT*,  STATUS  -  »OLD‘) 

READ  (1.SOO)  IFSR 

READ  (1,600)  (IFBO(K}  .  K  «  1,  IFSR) 

500  FORMAT(13) 

600  FORMAT (13) 

CLOSE  (UNIT  -  1) 

DO  S  K  *  1  ,  IFSR 
IVAL(K)  «  1 
5  CONTINUE 

DO  100  X  «  l.N 

C  CHECK  TO  SEE  IF  WE  NEED  TO  GENERATE  A  NEW  DATA  BIT 
IP  ( (K+TOELAY) *TSAHP  .LT.  I*TDATA*TSAMP)  GO  TO  10 
C  IP  SO  DO  IT  HERE 
RANDOM  «  RAND(L) 
c  debug 

C  imiTE(6,«)  *  IFLAG  AND  RANDOM  NUMBER  RETURNED: '.L, RANDOM 
INFO  •  0 

IF  (RAND(M  .GT.  0.5)  INFO  «  1 
C  KEEP  TRACK  OF  WHICH  DATA  BIT  HE  ARE  ON 
I  *  141 
10  CONTINUE 

C  HOW  CHECK  TO  SEE  IF  WE  HEED  TO  GENERATE  A  NEW  CHIP  BIT 
IF  ( (K4TDELAY) *TSANP  .LT.  J*TCHIP*TSAMP)  GO  TO  20 
IP  (ICHIP  .HE.  1)  GO  TO  15 
CHIP  -  IVAL(IFSR) 

CALL  MSEQ(IFSR.IFBD.IVAL) 
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GO  TO  19 

15 

CONTINUE 

C 

IF  SO  DO  IT  HERE 

RANDOM  *  RAND(L} 

CHIP  «  0 

IF  (RANDOM  .GT.  0.5) 

CHIP  *  1 

C 

KEEP  TRACK  OF  THE  CHIP 

BIT  NUMBER 

19 

CONTINUE 

J  -  Js  l 

20 

CONTINUE 

C  NOW  WE  DETERMINE  THE  PHASE  SHIFT 

PHASE  «  FLOAT (MOD (INFO  +  CHIP, 2)) 

C  COMPUTE  THE  ARGUMENT  FOR  THE  COSINE  FUNCTION 

C  YARG  -  2.0*PI*(1.0/(TCARR*TSAMP))*(K+TDELAY)  +  PI*PHASE 

YSAMP(K)»  MAG  *  COS (YARG) 

*  100  CONTINUE 

C  NOW  SAVE  THE  VALUE 

c  ranga  OPEN  (UNIT=l, FILE* * bpsk. dat ' ,STATUS=' NEW ) 

C  CALL  DATE(BDATE) 

C  CALL  TIME(CTIME) 

c-ranga  WRITE(1,125)N,MA6,'TCARR,TCHIP,TDATA 

125  FORMAT(I5,5X,F4.1,10H*OSBPSK+M2,F7.2,lHf,F7.2,lHc,F7.1, 

*  IHd) 

c-ranga  WRITE (I, 150)  (YSAMP(I) ,I»1,N) 

do  121  1*1, n 
c  writed,*)  ysamp(l) 

121  continue 


ISO  FORMAT(8F16.e) 
CLOSE  (UNIT-1) 
return 
END 


t 


V 


SUBROUTINE  KSEQ(IFSR, IFDB, IVAL) 

C  THIS  SUBROUTINE  PERFORMS  THE  SHIFTING  OPERATION  OF  AN  IFSR 
STAGE 

C  FEEDBACK  SHUT  REGISTER,  WITH  FEED  BACK  CONNECTIONS  AS 
INDICATED 

C  BY  IFDB.  IVAL  IS  THE  INITIAL  CONTENTS  OF  THE  FSR  AND  WILL 
CONTAIN 

C  THE  FINAL  CONTENTS  AFTER  SHIFTING. 

C  W.R.  TUCKER  4  AUG  83 

INTEGER  IFDB(IFSR),  IVAL(IFSR) 

ISUM  -  0 

DO  10  I  -  1,  IFSR 
ISUK  •  IFDB(l)  «  IVAL(I)  'I-  ISUM 
10  CONTINUE 

IBIT  -  MOD(ISUH,2) 

DO  20  I  -  1,  IFSR  »  1 

IVALdFSR  +  1  -  I)  -  IVAL(IFSR  -  I  ) 

20  CONTINUE 
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IVAL(l)  =  IBIT 

RETURN 

END 

c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


ee 
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